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ON THE LINEAR BIRTH AND DEATH PROCESSES 
OF BIOLOGY AS MARKOFF CHAINS 


ANTHONY F. BARTHOLOMAY* 


HARVARD UNIVERSITY, 
BOSTON, MASSACHUSETTS 


Stochastic Markoff models for the linear birth and death population 
growth processes of biology are constructed using the Q-matrix method of 
Doob. The relationship of the stochastic theory to the classical deter- 
ministic foundations of these processes is stressed by showing in detail 
how the classical postulates are mathematically transformed via the Q- 
matrix elements into the basis for a stationary Markoff process with con- 
tinuous time parameter and denumerably many ‘‘population states.”’ It is 
shown that the resulting stochastic models predict that the population 
size will fluctuate about the deterministic time curve, the extent of fluc- 
tuation being measured by the variance functions. General formulas cover- 
ing all possible transitions from one population size to another are derived, 


1. Introduction. In recent times the classical deterministic 
theories for biological population growth processes have been aug- 
mented by the construction of stochastic models which have fur- 
nished a means for interpreting random fluctuations in the popula- 
tion size with time. W. Feller (1939, 1950), D. G. Kendall (1948, 
1949), and others have in this way emphasized the role of chance 
events in the determination of the size of a population of cells or 
organisms. The linear ‘‘birth and death’’ processes, the subject 
of this paper, are examples of the basic stochastic processes 
which have evolved. They have also served as mathematical 
models for related phenomena in other fields. 

These stochastic growth models are usually constructed from 
axioms based on generalizations of those of the simple Poisson 
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process. It is the purpose of this paper to identify the mathemati- 
cal bases of such processes more completely with the more general 
theory of Markoff processes. This is accomplished by expressing 
the basic assumptions of the growth phenomena directly in terms of 
the Markoff theory and constructing stationary Markoff models for 
them. 

J. L. Doob’s (1953) work on Markoff processes has been found 
most convenient in this regard. Using his mathematical theory, the 
representation of the growth process as a stationary Markoff chain 
is begun by translating certain basic properties of the classical 
deterministic description of the phenomenon into the constant ele- 
ments of a Q-matrix. The Q-matrix forms the coefficient matrix of 
two matrix differential equations, the ‘‘forward’’ and ‘‘backward’’ 
equations, which together with appropriate initial conditions deter- 
mine the basic matrix function P(t) composed of the homogeneous 
transition probabilities of the chain. Once the Markoff model is 
specified in this way, the g’s in particular applications may be 
estimated from empirical data, using the maximum likelihood method 
or a modified y? method in statistics. Such statistical considera- 
tions have appeared in the literature (see Grenander, 1950; Lange, 
1956) which is almost devoid of estimation theory for stochastic 
processes in general. It is therefore useful to identify the growth 
processes with the only chapter of statistical estimation of random 
functions which has received much attention. 

The derivation of these processes as extensions of the Poisson 
theory is in some ways simpler, but the probability functions ob- 
tained are likely to be narrowly interpreted as the first-order dis- 
tributions of the process, when actually they are identifiable with 
transition probabilities in the Markoff sense, which, in turn, com- 
pletely determine the whole random function (that is, all higher= 
order probability distributions are obtainable from them). Incorpo- 
rating the Markoff features into the basic axiomatic construct of 
such processes is therefore important for both practical and theo- 
retical reasons. 

Thus, in the present paper, it will be assumed from the outset 
that the linear birth and death processes are stationary Markoff 
chains {v,, 0 <¢<} (see section 2) with continuous time param- 
eter ¢. The random variable z, has at any time ¢ the possible val- 
ues 0, 1, 2,..., which completely define the basic ‘‘population 
states’’ Sy, S,,..., into which and out of which the population 
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moves in accordance with the transition probabilities Pit) from 
state S, to S;. These ideas are further explicated in the following 
section where the pertinent basic concepts and theorems of the 
Markoff theory, following Doob (1953), are given. 


2. The theory of stationary Markoff chains with continuous time 
parameter and denumerably many states. The probability of pas- 
sage of a system from state S, to state S, t=O, 1 os. es) in the 
time interval from ¢, to ¢, + ¢ is independent of ¢, and is given by 
the ‘‘transition probability function’”’ p,,(¢) defined as the condi- 
tional probability Pla, = S; | t,. = S,} (for all ¢,) that the random 
variable z, has the value S,, given ¢ units previously it had the 
value S;. For the process {z,, 0 <¢ <0} to qualify as a stationary 
Markoff chain the transition probabilities are required to satisfy 
the following conditions: 


Pp, (2) 2 0 (i, j= 0, 1, Dsleie es) (1) 
>, Rift) = 1 (¢ = 0, 1, 2,...), (2) 
j 
(0<s<t<o), 

ye D,;(8)° P jx (4) = P48 + t) (3) 

j (PE 1 Poe FSS). 

It is also assumed that 
‘ sue Lan (te, 4) 4 
Lin f= 84 = 9 Ht 2 
and that a matrix Q = (q;;) of constants is defined by the limits: 
; .(t) -— 1 ‘ 
q;; = lim a = pj; (9), 


‘ P ijt) , 
sien Pat PLE 


These conditions and assumptions lead to the following systems of 
differential equations,* in which pi ;(2) stands for the time deriva- 
tive of P, jd): 


FORWARD SYSTEM p’(4) = du, Pin +), Qn Pi — (8) 
i#k 


*See Appendix. 
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BACKWARD SYSTEM 7%, (2) = q,:P;q (4) +). 9:3 Pix (2 (7) 

jAi 
In determining the transition probabilities from theoretical con: 
siderations: the q;; are first determined, g,;,; is set equal to 


- ay Qi; » an initial condition is postulated, and the resulting 
j#i 
systems (6) and (7) of differential equations are solved. In case 
the number of states, N +1, say, is finite, with the initial condi- 
tions Pp, (0) =1 (i=), P;;(0) = 0 (j #7), it may be shown (see 
Doob, 1953, pp. 240-41) that a unique solution to both systems is 
obtained in convenient matrix form 


P(d) = ef", (8) 


where Q is the matrix (q;;) and P(t) is the matrix of all transition 
probabilities [p,,-(¢)] and e**® is the element-by-element sum of 
the exponential series expansion: 1 + ¢Q + t7Q7/2!+.... 


In the more general case, where the number of possible states 
is infinitely denumerable, the situation is more complicated, and 
the discontinuities of the ‘‘sample functions’’ must be studied to 
determine whether a unique solution to both systems does exist. 
This rather involved question has been studied by Doob (1945, 
1953), and by Feller (1940, 1950) and will not be discussed here. 
Suffice it to say, if q;; (, 7,20; @# 7) is a sequence of non- 


negative numbers with g;; = — ye q;; for all 7, then it can be shown 
T#i 
that there exists a corresponding stationary Markoff transition 
matrix function P(t) = [p,,(¢)] which satisfies the backward sys- 
tem (7); the forward system (6) must be altered by changing the 
equality to the sign >. A modified matrix function P(t) = (7, ,(91 
in which only transitions S$; —> S, in finitely many steps are al- 
lowed can be found, however, which does satisfy both (6) and (7), 
Subject to the conditions 7; (0) = Pp; (0) =9,,, and p,(t) < P; (4. 


In the event that the series » rE diverges, only transitions 
i 

S; — S; in finitely many steps can occur, in which case P(2) = 

[p, (2 would be an unrestricted common solution. Thus the diver- 

gence of the reciprocal q;; series will in general guarantee the ex- 
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istence of a common solution to both (6) and (7). This will be 
seen to be true in the case of the birth and death processes de- 
rived here. 


3. The linear birth process. Assuming a population of organisms 
which increases monotonically in size as the result of binary fis- 
Sions, say, and with no deaths occurring, its rate of growth, ac- 
cording to the *‘Law of Malthus,”’ is given by the equation 


= =Az, (9) 


where A is a constant, and z the number of individuals present at 
time ¢. Where 2, is the number present at time ¢ = 0, integration of 
(9) gives 


a(t) = ayer. (10) 


According to the deterministic theory, this expression is sufficient 
to predict perfectly the population size at any time, once the ini- 
tial size and the growth constant A are specified. 

Actually, the size of such a population ranges over the set of 
integers, so that the continuous function representation 2(¢) is 
merely a mathematical convenience. Furthermore, in view of ob- 
served irregular fluctuations about such values, it seems reason- 
able to treat the population size as a random time variable and to 
propose as a simple stochastic model which would account for such 
fluctuations a stationary Markoff chain {z,, 0 <¢<e}, in which 
the ‘‘states’’ of the linear birth process S$,,S,, ..., are defined as 
the population sizes 1, 2,..., respectively. 

The specification of the stochastic process is suggested by the 
following information contained in the deterministic result (10). In 
the time interval (¢,, ¢, + Az,), the ratio of the number of births 
Az, to the number 2(¢,) = x, of individuals present at the begin- 
ning of the interval is given approximately by 


Az, 


vy 


= At, + 0(At,), (11) 


where 2, = a ers and o(AZ,) is the sum of infinitesimals of higher 
order. This expression may be interpreted as the probability of an 
individual binary fission in an arbitrary time interval of length A¢,. 
Assuming then that each of the z, present at ¢, has this same 
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chance of reproducing, the probability of exactly one such fission 
from these 2, individuals would be [Az,Az,+ o(AZ,)]. This fol- 


Ly cal 


lows by taking p = and g=1- in the Bernoulli formula 


b(1; z,, p) for 1 success in 2, trials. By analogy, in the language 
of stationary Markoff chains, let it be assumed that in a time in- 
terval of length AZ, (independently of the beginning point) only the 
transition S, + S,,, out of state S; is permitted, where now 7 is 
the population size at the beginning of the interval. That these 
considerations determine the Q-matrix may be seen as follows. 

Expanding the transition probability function DP; ;(o) about 0 as a 
Taylor series [assuming that p,,(¢) is analytical] gives 


Ad? 
Pp, 42) = p, (0) + p,(0) Az + p7i(0) perens ue (12) 


Thus, setting 7 = 7+ 1, one obtains 


Pigia1\A2) = G3 541° 40+ O(Ad) (13) 
which by definition of the transition probability function is the 
probability of a transition S; — S,,, over a time interval of length 
At. However, according to the previous paragraph, p, i41(A t) 
should have the value [AiA¢ + o(Ad)]. Thus, 


Qi, i41 = Xt @=1, y. oer pe (14) 
In addition, let 
qi; = 0 (j4#7%,2+1), (15) 
and 
Gi Via Nt (§=1,2,...). (16) 


This completes the specification of the Q-matrix and leads to the 
systems of forward and backward equations by substitution of the 
appropriate q-values into equations (6) and (7): 
FORWARD SYSTEM 
OF THE LINEAR P(t) = —kdp,, (t) + (K-1)A 

k Vik P; — (2) (17) 
BIRTH PROCESS: nape” 
BACKWARD SYSTEM 


OF THE LINEAR 9’, (t) = -1Ap,.(#) + iAp,., (2). 18 
BIRTH PROCESS: s ® ttn) Sipe 
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é = ary et te RES 

Since ys Goo = = De (~') diverges, a unique solution P;,(¢) of 
t=1 t=1 

both (17) and (18) is known to exist. This solution may be ob- 

tained from (17) by induction because of the recursive nature of 

this system of equations. A more generally applicable method 

which makes use of the ‘‘probability generating function,’’ 


oo 


peal. (p,,(2)l = $;(8,2) =D) 3* p;,(2) (19) 
k=1 
(where s is an arbitrary variable), will be used to solve (17). 


Differentiating ¢,(s,¢) partially with respect to ¢ and substituting 
from (17) into the result gives 


ad, = ~ 
oF Piss k : 
Soave aise oe ks ‘pelt +>) (k-1)s Ps ~—164)3 (20) 
k=1 k=1 
99; ==As E ks*~*p,,(t) +2 e ks**" p(t), (21) 
at ; 
k=1 k=1 
or 
ab, _ eee 25) pak! 29 
oie -—r8 » ks*~p,,(t)+r8 ks*~" p40) (22) 
: k=1 kek 
A) . co 
But a =) ka" (23) 


so that (22) reduces to 


0 


“St =n o(s-1) 


a¢, 
2 (24) 


This is a linear, homogeneous, first-order partial differential equa- 
tion, the solution of which may be obtained using the auxiliary 
ordinary differential equation method, in which (24) leads to the 
auxiliary equation: 


(25) 
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the solution of which may be expressed as 


s-1 
Ga eee eee (26) 


s 


where a is an integration constant. Thus, the general solution of 
(24) may be expressed as 


$(s,t) = O(a) = 0 (- ? : 2%) (27) 


The fact that S, occupies in p;,(¢) the role of initial state becomes, 
in terms of ¢,(s,¢), ri Z 


$,(8,0) = si= ® (- : *) (28) 


and provides a boundary condition for the partial differential equa- 
tion (24). A second boundary condition which will be automati- 
cally satisfied due to the nature of the Q-matrix is 


$(1,)= ) > pO =1 (29) 


k=1 


which is obtained from (2) in section 1. 
Introducing the new variable 


s-1 


s 


S= 


, (30) 


allows ® from (28) to be expressed as 
G(s) = (t=s)-', (31) 
from which it follows that 


s-i 


,(s,t) = ® ( et) = [1 —(1 = 74) s)78gte- i, (32) 


Expanding the right-hand side of (32) gives for the (k — 7+ 1)th term 
(# me ‘) (ars e7 M)k-i ew dit gk. (33) 
—i 
Comparing this with (19) leads to 


P;, (0) = o : :) (i - e Mtyk~i ew hit, (34) 
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where 


a k-1 ' 
(1 <i<k<o), or p,,(t) = bay en Mkt (prt _ 1)hmi, (34) 


As a special case of (34) or (34'), 


P;;(¢) = ens (35) 
indicating that the probability of no transition out of state S, de- 


cays exponentially, passing to the limit 0 as ¢-+ ~. In fact, it is 
generally true that 


lim p,,(t)=0 (k> 4), (36) 
t 20 


indicating that there is, after a long time, little chance of a ‘‘small”’ 
or finite population. On the other hand, 


lim p, (2) = p;{0) = 1, 
t>0 (37) 
Pinft) = Py, (0) = 9 (A> 2). 


The mean value function and the variance functions of this proc- 
ess are most easily calculated from the formulas: 


Od. 
= ), (38) 
s=1 


d8 
+ E(t) - EX(d). (39) 


Mean Value = E ,(¢) = ( 


a? 
Variance = o?(t) = (24 
s 


Since from (32) the required partial derivatives evaluated at s =1 
have the values 

Od; 

08 


2 
0d, 
ds? 
it follows that the mean value, corresponding to S; as initial swte, 
is 


= ier, (40) 


s 


wie [(i + 1)e** - 2], (41) 


s=1 


E(t) =ie™, (42) 


and the variance 
o3(t) = ie™ (e™ - 1). (43) 
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In the notation of the deterministic equations, (9) and (10), if 
the initial state S; is identified with a, and the state S, with z 
then the probability of z at time ¢, from (34) becomes 


P(e) = G 7 Ja e7 *Xot (1 - ew M)x—xo, (44) 


The mean value (40) in this language becomes 


E(t) = z,e™, (45) 
and the variance (43) 
o2(t) = xe (e™ — 1). (46) 
But the right-hand side of (45) is the same as the right-hand side 
of (10). Thus the foregoing stationary Markoff process describes a 
fluctuation about the deterministically predicted values, a measure 
of the dispersion being given by the square root of the expression 
(46) for the variance. Furthermore, equations (44)-(46) are the 
same as those obtained for the corresponding Poisson-type proc- 
ess (see Kendall, 1949; Feller, 1950). 

That (84) is a solution to both the forward and backward equa- 
tions, (17) and (18), respectively, may be verified by direct substi- 
tion. Thus, the identification of the simple linear birth process as 
a stationary Markoff process is complete. Because of this identifi- 
cation the stochastic process is completely specified; viz., all 
higher-order distributions are determinable from the transition 
probability matrix function P(t) = [p,,(z)]. 


4. The linear birth and death process. Suppose that a linear 
death force, i.e., one proportional to the population size, is added 
to the assumptions of the linear birth process of the previous sec- 
tion. In classical deterministic terms this is dealt with by setting 


dz 
the rate of growth F of the population size x equal to the differ- 


ence between the rate of increase due to the binary fission process 

and the rate of decrease y « due to the deaths of individual organisms: 
dz 

a 7 hee =(C-a)e. (47) 


Together with the initial condition 2, = 2(0) this completely de- 
termines the population size at any time: 


t= 2 e(h— Ht, (48) 
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Thus, the deterministic mathematical treatment of the birth and 
death process is no different from that of the simple birth process, 
the difference 0 between the birth and death constants, 0 = A - My 
replacing the simple birth constant d. 

The problem of constructing a stochastic model for this process, 
which is capable of providing a mathematical description of the 
fluctuations about the deterministic curve (48), is more involved 
than in the case of the simple birth process, however. The birth 
and the death components must be handled independently, not 
simply as a difference. 

Making the assumption that the linear birth and death process is 
a stationary Markoff chain {z,, 0 <¢ <0}, a Q-matrix may be con- 
structed from the following inferences derived from the determinis- 
tic. model. From (48) it follows that the ratio of the net change A a, 
in the population size to z, = #(¢,) over the interval (¢,, ¢, + A¢,) 
is given by 


Az 


x 


~=AAt, -pAt, + 0(At,). (49) 
1 

But where Ad is the number of binary fissions taking place in this 
time, and Ad the number of deaths, it follows that 


=-—-—. (50) 


The first term on the right gives approximately the probability that 
an individual subdivides by binary fission in this time; the second, 
the probability than an individual organism will die. 

The probability of exactly one birth or of exactly one death from 
the z, individuals present at time ¢, in the time (¢,, ¢, + At,) now 
may be easily calculated by assuming that %, ‘‘trials’’ are made, 
each resulting in one of three possible outcomes: (1) reproduction 
by the individual, (2) death of the individual, (3) neither a reproduc- 
tion nor a death. Then the multinomial distribution law with p, = 
aed =)At, + 0(A2), po= ae =pAt, + o(At,), Pg =1- (p, + Pa) 
ick: that the probability éf exactly one birth in this time is 
[Aw At, + o(AZ,)], and the probability of one death equals [ua,At, + 


o(At,)). 
Nits that in this analysis the probability of more than one birth 
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or more than one death receives the asymptotically negligible 
value o(A¢,). In other words in this time the population is not 
allowed to pass from size a, to (2, + 1) by a difference of 1 be- 
tween an arbitrary number of births and deaths. Incorporating this 
feature into the stochastic model then leads to a very simple model 
compared to what might be obtained from probability assumptions 
permitting multiple events in the interval of length Az,. The fact 
that A¢, is an infinitesimal in these considerations makes the 
simpler assumption a reasonable one. 

The population states S,, S,,..., of the Markoff chain will be 
defined as in section 3. In addition, the state aes corresponding 
to extinction is added and is referred to as an ‘‘absorbing state,”’ 
connoting that once the ‘‘system’’ reaches this state it remains 
there. Making the same assumptions which were indicated in the 
previous paragraph, it is seen that in an arbitrary time interval of 
length Az (independent of the beginning point) only the transitions 
8S, — S,;_,, 8; + 8;,,;(@=1, 2,...) are permitted and that the 
transition S, —> §, receives a zero probability. Thus from (12) it 
follows that 


P;,:-1(A%) = 9;,;_, At + o(Ad), (51) 
Ps ig (AD = V:,ig, 40+ o(At), (52) 
and hence that the Q-matrix is defined by 
Tp jan. = ths 
qi,ia1 = OAs (fe O° 175s (53) 


EF -t(A +p), 


The divergence of the series >» (q;,)" ' guarantees a common solu- 
i=0 

tion to the following systems of equations obtained by substitu- 

tions of (53) into (6) and (7): 


Pin(t) = —k(A + 0) pi, (t) + (R- LAD; 4_ (+ 
(e+ 1m; 44 (Oy (54) 
Diy (t) = -2(A + n) p;,(0) + tADdi 41 4(D + TP: 1 p(s (55) 
(i ka0 eS 


respectively, the forward and the backward systems of differential 
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equations of the linear birth and death process. Where (1 <i < k < «) 
and » = 0 equations (54) and (55) reduce to equations (17) and (18) 
of the simple linear birth process. These systems of equations are 
not in recursive form, so that the inductive method cannot be used. 
In this case, then, one is forced to use a method which determines 
the [p,,(2)] all at once. Such a method is that of the probability 
generating function, already used in section 3. 
Defining 
p.g.f. [p,,(2)] 2 ¢ (8,¢) zd ae s*p,,(2), (56) 
k=0 
one obtains 


ad, = = 
a =-(A+n)s es bat "ps {) +s? Ds (& - 1) s*-?p, ,_ (4) = 
k=0 k=0 
k=0 
Because 
Od; = 
sete >= ksta 1p, (0, (58) 


equation (57) may be rewritten in the form 


ha Pe siuyeetes 59 
Peay 1)(A8 - 2) appa. (59) 


the fundamental partial differential equation of the probability 
generating function for the linear birth and death process. Again, 
this is a linear, homogeneous, first-order equation, with auxiliary 


differential equation: 


Te ee (60) 
-1 (s--1)(A8 -p) 

the solution of which is obtained in the form 
s-1 

AS-— UL 


Casel. A#u; a= eh Mt. (61) 
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1 
Case 2. A=yp3 @-= pend: (62) 
8- 


where a and a’-are integration constants. 


Discussion of Case 1. The general solution of (59) with A#u 
may be expressed as 


$(8,t) = ®(a) = ® la conn], (63) 
Again, using the boundary condition 
s(oo- 0 @[ ST, Sie! 
one obtains, after the transformation 
s- (65) 
the form 
(9) -[ = | (68) 
Consequently, 


* eset s(homne ee | (67 
he-u X(e- BOs ny te | 


[ue A-M* — yu] — [py eA-# a] 8]! 
soo [Reocmaactnomrcag) 


Expansion of the right-hand side of (68) in powers of s leads toa 
coefficient for s* which when matched with (56) gives the following 
values for the transition function p,,(t): 


(Hiak yn sont fay eat be fe 
Pik eet i n hmm, i . 


(69) 
[eCA-me L]itk-2n[y e(A-M)t pl-i kta Ly eC Ame _ Al", 


where the upper summation limit is & if 0 <k <i; and 2, if R27. 
As a special case of (69) when pz = 0 one obtains formula (34) for 
the pure linear birth process. From (69) it follows that 


Lim Paull) = Pal) = 5 EO (70) 


BIRTH AND DEATH PROCESSES AS MARKOFF CHAINS 111 


indicating satisfaction of the Markoff relations (4). 

The mean value and the variance functions are most conveniently 
calculated from the probability generating function (68), using for- 
mulas (38) and (89). Since 


ad, mee 
(= ) ri teas (71) 
sS=1 
and 
07d. z 
as He! (A=W trey _ (A=) ¢t 
(es =A oot. [e(A u)e ts 
, (72) 
(A + u) e' —M)t _ 2 dl, 
the mean value 
E(t) = te“ (73) 
and the variance 
r 
o*(t) =¢ — e(A-Mt [e(A-wet _ 4], (74) 


In the language of the deterministic theory [see equations (47) 
and (48)! if the initial state S, is defined by x, individuals, and 
the state S, by 2 individuals, then the probability of 2 at time ¢, 
call it p(t), is written: 


i acl x %+2,-n—-—1 
fas _4\% 0 0 Ax X0-n , 
p (2) D (-1) Gaal a ) mn 
r= 


(75) 
[eCA-Met _ 4] *0+2-2/y eA-Met _ p)-*0o-* +" [y eOh-Mt _ Al, 


and the mean value becomes 
A=) 
Eetese (76) 


Comparison of (76) with (48) shows the right-hand sides to be 
identical. Thus, the stochastic linear birth and death process 
gives a mathematical representation of the fluctuations of a popu- 
lation about the deterministically predicted values, with variance 
equal to 


r 
7 3(2) See : +p eA-HEf eA ws sii (77) 
—U 
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As explained at the beginning of section 4, a transition into 
state S, represents extinction of the population. The probability 
of extinction starting from any state S; may therefore be obtained 
from formula (69) by setting & = 0 and rewriting it in the form: 


A- pu : 
P iol?) = ; eh | : (78) 
1 (A <n), : 
s = 9 
It follows that lim P ;oft) 0 (sw. (79) 


so that, if \ <u as a limiting case (see also Case 2 following), 
extinction is ultimately certain, in a probabilistic sense. This is 
a result that could not be predicted in the deterministic theory. 
Discussion of Case 2. X=. The ideas of the preceding dis- 
cussion of Case 1 (A #4) as given in formulas (63)-(79) may be 
applied to the particular case A = yu (i.e., in terms of the determin- 
istic theory, birth rate constant equals the death rate constant). 
As already shown, the corresponding partial differential equation 


ALINE eet (80) 
Ot ds 
leads to the ordinary auxiliary differential equation 
dt ds 
ae Neb? (81) 
with solution expressible in the form 
1 
@’-= -Xt, (82) 


8-1 


a’ being an integration constant. 
Thus, in this case the general solution of (80) is expressible as 


1 
blot) = (a) = 0( > -at). (83) 
And, using the boundary condition 
: 1 
¢ (3,0) =s'=9 (|, (84) 
: 1 
the transformation S = (85) 
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Saci\* 
so that 
: hea _ faAt-(Qe-1)s]! 
OS heey (- =| at) ‘ E Nd) 2] (87) 


gives the probability generating function, the solution of (80). Ex- 
pansion of (87) in powers of s and in comparison with 


 ,( 8,2) = 3 s* p.,(t) 
k=0 


then gives the transition probabilities for the (equal birth rate- 
equal death rate) birth and death process: 


Pix(t) = as cayr(§) ( begs 82 ; 
n=0 


(At) tt—32 (At = 1)” (At . ieee 


(88) 


the common solution of both the forward and backward systems of 
differential equations: 


F.E. pi, () =-2hAp,,() + (A - 1A, ,_,() + (39) 
(k + WAP; pa 1s 
B.E.  pi,(t) =-2¢Ap,,(t) + trp, (2) + tADs_ 1, (0- (90) 


Equations (87) and (88) may be obtained from (68) and (69) by pas- 
sage to the limit as » —» ». Independently of Case 1, from (88) it 
may be seen that 
Pio(@) = (Ad)(1 + At, (91) 
from which it follows that 
lim p,,(¢) = 1, (92) 
t 00 


so that extinction is ultimately ‘‘certain’’ in this case. It is also 
true in this case that 


‘3 (é =k), 
lim P ix (t) =0 (i # ah : (93) 
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Furthermore, since from (87) (0¢,/0 8220 =, (94) 
and (P4/ds7| _,)=t(2rAt+2- 1), (95) 
it follows that 
E(t) =t, (96) 
o2(t) = 2iAt, (97) 


respectively the mean value and variance functions of the special 
process. As (96) and (97) show, in this case as the growth proc- 
cess continues there is a fluctuation about the initial size ¢ which 
increases indefinitely with time. 


5. Conclusions. The classical linear birth and death processes 
of biology have been re-examined strictly from the point of view of 
the mathematical theory of Markoff processes, using Doob’s Q- 
matrix method of synthesis. It has been shown that these biologi- 
cal processes may be represented by stochastic models which are 
stationary Markoff chains with continuous time parameter and de- 
numerably many states. The probability of a given number of indi- 
viduals at any time is interpreted as the probability of transition 
from some other state to that state occurring in a time interval 
equal in length to that time. These transition probabilities [p,,(2)] 
define completely the random growth process. 

Many of the results of this investigation appear scattered in 
various sources listed in the references. In this paper a more com- 
plete discussion has been given and the unifying factor has been 
the straight Markoff approach. The results obtained by Feller and 
Kendall using the Poisson method are shown to be a special case 
of the formulas derived. The birth and death equations appear in 
the literature only for the case in which the process begins with a 
single individual (corresponding to t=1). In this paper this re- 
striction had been removed and the results obtained therefore ap- 
ply to any initial state S,, where z > 1. 

Particular attention has been given to the problem of stochastic 
model building by giving in detail the mode of transition from the 
deterministic to the stochastic theory. The birth and death rate 
constants are carefully traced from the deterministic context to 
their role in the determination of Doob’s Q-matrix. This matrix of 
constants which comprises the coefficients of the backward and 
forward differential equation systems is the transformer from the 
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deterministic to the stochastic model, for once these systems are 
solved the transition probabilities [p,,(¢)] and all higher-order 
distribution functions of the process are automatically determined. 
The fact that the resulting probability functions have mean values 
which are identical with the deterministic equations shows that 
the particular mode of construction of these stochastic models has 
led to models which describe fluctuations about the values ex- 
pected according to the earlier theories. The corresponding vari- 
ance functions then give an idea of the range of these fluctuations. 

The question of ultimate extinction in the case of the birth and 
death process is one which cannot be discussed in terms of the de- 
terministic model. In the case of the Markoff model, ‘however, the 
problem of ultimate extinction is discussed in terms of transitions, 
S, — S,, into the ‘‘absorbing state’’ consisting of 0 individuals. 
It has, in fact, been shown in section 4 that with increasing time 
extinction grows asymptotically certain when the death rate con- 
stant 1 exceeds, or equals, the birth rate constant; when the oppo- 
site is true, the probability of ultimate extinction has been found 
to be the ratio of the death constant to the birth constant raised to 
the power of the initial size. 

In the case of the pure birth process, discussed in section 3, it 
has been shown that lim p;,(t) = 0, indicating that in sufficiently 

00 


long time finite states receive a low probability. From ergodic 
properties of the Markoff chain theory one would also expect this 
property to be true in the case of the birth and death process. In 
this event one could then describe the process by stating that after 
a long time the population, starting from any (non-zero) initial state 
becomes extinct, or else tends to infinity, the probability of transi- 
tions into finite states approaching 0. Using the general formula 
(69) we have been able to check for various values of z that this is 
true, but the general proof has not been obtained. 

A difficulty which arises in representing phenomena of this type 
as stochastic processes is the solution of the systems of equations 
which are produced. The method of the probability generating func- 
tion, advocated by Feller in the treatment of such problems, allows 
one to transform such systems into partial differential equations. 
In the case of the birth and death processes the partial differential 
equations are first-order, homogeneous equations which yielded to 
the method of the auxiliary differential equation. The main details 
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of this very useful method have been given, and it is hoped that 
they may be of some use to others faced with similar problems. 

While the direct application of these processes are to the field 
of biology, it should be emphasized that such processes have a 
very wide range of application. Modifications of them have, for 
example, been employed in the study of fluctuations in cosmic 
radiation (see Arley, 1943), in the queueing problems in telephone 
communication theory (see Feller, 1950), and quite recencly by the 
author in the study of chemical reaction kinetics. 


I should like to thank the editors for their helpful suggestions, 
all of which were followed in the revision of the original manu- 
script. I owe a particular debt of gratitude to Dr. H. G. Landau 
who called my attention to the results obtained from Equation (78) 
which had been incorrectly stated by me in the revised copy, and 
who made other useful comments. 


APPENDIX 
Derivation of equations (6) and (7) of the text. 
1. Starting with Equation (3) of the text: 
Pix(8 + t)= Dp, (8)+D;(2) 
j 
(where 0 <8<t<oo), the ‘‘forward’’ system of equations is ob- 
tained as follows: 
ee +t) 
ik 22 P (8+ Pi (t)- (1A) 


In this and the following equation p’ AAC) stands for differentiation 
with respect to time. Thus, 


Op, (8 + ¢) 
ee ieee OM Ae (2A) 
dp. 
But P ,(8 + t) taapes Dix(8 + At) —p,,(8 + tbe (3A) 
Ot t=0 At 0 At 
ey Pi,(t + At) —p,,(t) (4A) 


At>0 pare 
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i.e., 


Op .,(8 + t) 


ay = D;,(t). (5A) 


t=0 


Rewriting the right-hand side of (2A) and introducing (5A), the re- 
sult obtained is the ‘‘forward system,’’ Equation (6) which we re- 
peat here: 


Pil) = Ven Pin + D> 9), P12 (8) 
j#k 


2. Differentiating (3) partially with respect to s gives 


Op,,(8 + t) Op ;.(s) 
SE OIC eer ee (6A) 
j 
whence 
Op ,,(8 + t) 
ds = x 2 P ;,(¢) Vij (TA) 
But 
PP(e+O) —_  PiglA s+ t)- 7,,(0 + 0) (8A) 
ds jug 80 As 
ik(t+ At)-p.,(t 
mp paBdL VaR nUe (9A) 
Ato At 
or 
OD. t 
Pix(8 ue ) = D4 (t)« (10A) 
0s 23 


Substituting (10A) into (7A) and rewriting the right-hand side, 
gives then the ‘‘backward system,’’ Equation (7) of the text: 


(AC) = 9:;P ix a ee 9; ;P j4()- (7) 


j#i 
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This article gives a general mathematical description of the develop- 
ment of a colony of social wasps. It also gives a mathematical deriva- 
tion of the time when it is most suitable for a colony of wasps to begin 
to produce queens (instead of only workers). First, we consider the case 
when the wasps do not build special cells for the queen-eggs and then 
the case when they do. By comparison with observed data which are 
available in the literature the theory seems to give a very good descrip- 
tion of what actually happens in a colony of social wasps. The theory 
can be used to check the validity of certain hypotheses trying to explain 
when and why a colony of wasps suddenly begins to produce queens. 


Introduction. As is well known the males and the workers among 
the social wasps in most parts of Europe and North America, or 
generally on the parts of the earth which have cold winters, usually 
die in the fall and only the young fertilized queens hibernate. In 
the spring the queens wake up from their winter doze and immedi- 
ately begin to build their nests. The eggs a queen will lay in the 
beginning will all develop into workers. Presently the colony will 
grow, the workers will build more cells and combs, but still there 
will only come workers from the eggs. Suddenly, however, the 
colony will begin to produce males and queens. Probably there 
will develop many workers now too but the primary importance is 
attached to producing reproductive individuals. Now a very na- 
tural question is: When do the wasp colonies begin this production 
of males and queens? It seems to be evident that there must bean 
optimal time to begin. If the wasps begin early, they will certainly 
have more time to rear queens but there will not be so many workers 
to feed and nurse the eggs and larvae; if they begin late there will 
be many workers but shorter time to spend. The purpose of this 
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article is to derive a mathematical theory for the development of a 
wasp colony and the optimal time which we can expect from mathe- 
matical considerations. By lack of accurate observations of this 
time we cannot compare the results with experimental values in 
detail but the results obtained here seem to be very reasonable 
and conform to the estimations given in the literature. 

As mentioned above, we will first treat the case when the wasps 
use the same sort of cells for the reproductive individuals as for 
the workers. This happens, for example, among many species of 
Polybia. We make the following assumptions: 

Assumption 1. When the production of reproductive individuals has 
begun, the ratio of the number of eggs laid per unit time which 
develop into workers to the number of eggs laid per unit time which 
develop into queens is constant. This can perhaps at first sight 
seem to be a startling assumption but we think it is justified for 
three reasons. The first, and most important, is that nobody seems 
to have measured this quantity and so we do not know how it 
varies. Therefore it seems rather senseless to assume some func- 
tional form of it. The second reason is that even if this ratio 
varies we can take an average value of it and use this value in- 
stead of the real value. This is what is often done in physics and 
other fields which apply mathematical methods. The third reason 
is that anyone who wants to use some special functional form of 
the ratio other than a constant can do that by only a slight modifi- 
cation of the treatment in this article. 

Assumption 2. In the second case, when the wasps build special 
cells for the queens (species of Vespula, for example) the main as- 
sumption is that the ratio of the number of cells built for queens 
per unit time to the number of cells built for workers per unit time 
is constant. This assumption can be defended by the same rea- 
sons as in Assumption 1. 

The following assumptions are made in both cases. 

Assumption 3. The proportion of the (worker) cells which are 
empty at any moment is constant. This has been shown to be ap- 
proximately true for Polybia and the value is about 1/5 (Richards, 
1958, p. 71). 

Assumption 4. The number of cells one worker can build per unit 
time is constant. 

Assumption 5. The number of ‘‘nurses’’ which each egg, larva, 
and pupa requires is constant. Now of course a larva needs more 
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care than a pupa, but if we take an average value of the numbers 
which are needed for each of these three stages we obtain a value 
which certainly will give as good results as if we had assumed 
different numbers of ‘‘nurses’’ for each stage. 

The mathematical treatment given below is divided into three 
parts for both of the above cases. The first part deals with the 
state of the colony before the time when eggs which will develop 
into queens are laid. The second part deals with the time after 
such eggs have been laid but before any queens have come out of 
their pupas, and the third part deals with the time after the first 
queen has appeared but before the colony breaks down in the fall. 

The mathematical treatment will begin when the first workers 
have appeared in the spring. It would be of no use to begin earlier 
because everything is too uncertain before there is a definite num- 
ber of workers in the nest. We must not forget that a theory like 
that given in this paper is based on statistical assumptions (as 
those assumptions already stated) which are certainly not valid if 
there is just one wasp in the nest! As the lone queen builds be- 
tween 10 and 20 cells (Richards, 1953, p. 47), we can begin the 
mathematical treatment when the first 10 to 20 workers have ap- 
peared. This time is chosen as zero in the following. The results 
obtained below are valid for any number of egg-laying queens in 
the nest. 


We will now give a few definitions: 


n(t) is the number of cells existing at time ¢. 
n,(¢) is the number of workers living at time ¢. 
n,(¢) is the number of workers which are building new cells. 
n,(¢) is the number of workers which are acting as ‘‘nurses.’’ In 
this number we must include the workers which are out 
searching for food, the workers which are guarding the 
nest, and all other workers which have something else to 
do other than building new cells. 
n, is the number of ‘‘nurses’’ each egg, larva, or pupa needs. 
n,(t) is the number of cells which are empty at time ¢. 
n,(t) is the number of brood at time ¢ which will develop into 
workers. 
n,(t) is the number of eggs which will develop into workers, laid 
before and up to time ¢. 
is the number of cells one worker can build per unit time. 
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N(t) is the number of brood living at time ¢ which will develop 
into queens. 
N,(¢) is the number of queens living at time ¢. 
N,(t) is the number of eggs which will develop into queens, laid 
before and up to the time ¢. 
is the time which elapses from when an egg is laid until it 
has gone through the larva and pupa stages and the ready 
worker appears. 
t, is the average time which a worker lives. 
, is the time it takes for an egg which will develop into a 
queen to go through all the egg, larva, and pupa stages. 
We always assume in the following, when nothing else is 
said, that ¢, 2¢,. 
T, is the time when the first egg which will develop into a 
queen is laid. 
T is the time in the fall when the colony will break down. For 
tropical species 7 can instead denote the swarming time 
of the young queens or a similar quantity. 


Now we begin the treatment of the case when the wasps build no 
special cells for the reproductive individuals. 
PART ] 


Case 1. 0<t<T,. The equations for the state of the colony 
will be (their explanation follows immediately): 


ng(t) = ant), (1) 
n(t) = n,(t) +, (¢), (2) 
n,(t)=n,(t) +2,(2), (3) 
n, (¢) 
n,(t) = inka ; (4) 
ns (t) = n(t) — ng (¢) + n, (é) +n, (¢-t, - ts), (5) 
n,(t)=n,s(¢-t,) -n,(t-2¢, -t,), (6) 
n*(t)=ngn,(t), (7) 
g dn(t) 
where n’ (t) = mecrme the derivative of n(t) with respect to ¢. 
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Equation (1) states that the proportion on of the n cells is 
empty. According to Assumption 3, & is constant. 

Equation (2) states that the number of cells equals the number 
of empty cells plus the number of worker brood. 

Equation (3) states that the number of workers equals the number 
of ‘‘builders’’ plus the number of ‘‘nurses.”’ 

Equation (4) states that there cannot be more brood than the 
number of ‘‘nurses’’ allows. 

Equation (5) states that the total number of eggs laid up to the 
time ¢ equals the number of brood at time ¢ plus the number of 
workers at time ¢ plus the total number of eggs which were laid 
before time ¢ — (¢, + ¢,). All the wasps originating from these eggs 
are dead at time /. 

Equation (6) states that the number of workers living at time ¢ 
equals the number of eggs laid between the times ¢~¢, and 
t—(t, +¢,). The eggs laid later are not full-grown wasps yet. If 
we want to take the egg eating and the infant mortality into ac- 
count we need only multiply the right-hand side of equation (6) by 
a quantity less than one. In order not to encumber the analysis 
with too many constants we have not included this constant. 

Equation (7) states that there cannot be more cells built than the 
number of building workers allows. 

We will now solve this system of equations. It turns out to be 
most convenient to express all quantities in n, and then solve the 
equation for this function. From (1) and (2) we obtain 


n,(t) =(1- a)n(Z). (8) 
From (4) and (8), 
n, (t) =(1- %)n,n(Z). (9) 
Equations (7), (9), and (3) give 
n, (t) = = n’(t) +(1~a)ngn(d). (10) 
8 
Equations (5) and (6) give 
ny (t) = n(t) —n, (t) +75 (¢-7,); 


or with (1), 
n,(t)—n,(t -—t,) =(1— %)n(t). (11) 
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Equations (11) and (10) give 
1 


1 
n,(t) =—-: 

16 
From (12) and (6) we can finally obtain 


ji [nj (t) — ni (¢ —t,)] + ngng(d) —Ngn,(t — ¢,). (12) 


4 


ni(t)—n,(¢-—¢,)+n,n,(1 — o)n, (2) - (13) 

(n, +1)n,(1-4)n,(¢-¢,)+,(1-4)n, (¢-¢, —¢,) =0 
Expression (13) is an equation in n,(¢). It is known that such an 
equation has a solution of the form 


n,(t) = Ae**, (14) 


where A and 2 are constants. The value of z is determined by 
equation (13) and A by the initial conditions. Substitute the ex- 
pression (14) for n, (¢) into equation (13). After dividing each term 
by Ae** we obtain 


1 


e2-—2e “t+ nmon,(1 - 4) —(ng + 1)n, (1 - aje * + 


n,(1- aye *41t*? 20, fe 
This is an equation for obtaining the value of z. We see at once 
that one root is z= 0. We will show that this equation under very 
general conditions always has a root which is not zero. 

Call the expression to the left of the equality sign in (15) f(z). 
Then we know that f(0)=0. Further we can see that f(~) =o, 
f(-—~) =o. Because f(x) is a continous function of 2, the last 
equations mean that there are finite z-values, both positive and 
negative, for which f(z) > 0. Now, as f(0) = 0, a sufficient condi- 
tion that equation (15) has a root different from zero is that f*(0) # 0 
(see Figure 1). 

For f’(x) we find the following expression: 


f(z) =1-e "14+ te “14 t,(n, + 1)n,(1-a)e "1 — 


(¢, + t,)n,(1- a) got ag 
We find 


f°(0)=t,ngng(1- a) -tyn,(1- a). (16) 
The condition f’(0) 4 0 is then equivalent to 
t, #t,n,- (17) 
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FIGURE 1 


Even if ¢, =¢,n, it is not certain that equation (15) has only a 
root z = 0. It can have several roots even in this case. Now, how- 
ever, 7, is usually less than one and ¢, is quite certain less than 
t,. But if this is true, (17) is true too. Therefore we have not in- 
vestigated what happens when /’(0) = 0. 


We have now found that the solution of (13) is 
Ns (t) = Ae**, 


where z is a non-zero root of equation (15). 
From the last equation and (11) we obtain 


A —xt 
t) = e**{1 =e"). 
n(t) = ——— e**( ) 
The number of cells at ¢ = 0 is then 
(OP Se 715, 
1-4 


If we take as our initial condition the fact that there were n, cells 
in the nest at time zero we find the value of A: 


no(1 — a) 


: 18 
Oy ake aa 


A= 
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When we know the functions n, (¢) and n(¢t) we can easily find 
all the other functions too. The results for Case 1 are listed 
below. 


n(t)=n,e*', 


x 
nm, (t) = EB +(1- 2)n,| No age 
Ae 

2 
n, (t) = — 1, e*’, 

cs 


n (t) = N,N, (1 — a)e**s 

[ed 04 Sse geal el 

n,(t) =n, (l—aje**, 

n, (t) = Ae** A is defined in (18), 
N(t)=N, (0) =N, (4) = 0. 


From these expressions we see that the colony (the number of 
cells and workers) varies exponentially. Actually it will grow only 
if ze > 0 because otherwise the exponentials will decrease towards 
zero for increasing ¢-values. But ifn, <1 and¢, <¢, (as itis in 
nature) we see from (16) that f’(0) < 0, which means just that there 
is one value of e which is a root of equation (15) and which is 
larger than 0 (see Figure 1). 

It seems very probable that for each set of constants (&, ny, ng, 
me, ¢,, t,) there is one and only one way in which the colony can 
grow and prosper. That would mean that equation (15) has one and 
only one root different from zero. 

Now we can also see what would happen if condition (17) were 
not satisfied, that is, © 


to =t,Nn,- 


This relation means that the birthrate and the deathrate would be 
the same and that the number of workers and cells would be con- 
stant. If we inspect the results listed above we find that in this 
case «=0. That means that, if our mathematical description is 
correct, equation (15) will not have any non-zero root if t, =t, No, 
but will only touch the z-axis in the origin. 
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Case 2. T,<t<T,+t,. This is the interval in the life of the 
colony when eggs, which will develop into queens, have been laid 
but no such individuals have come out of their pupas yet. 

The equations for this case are 


n, (t) = an(Z), (19) 
n(t) =n, (t) +n, (t)+N(t), (20) 
n,(t)=n,(t) +n, (0), (21) 
n, (t) = 3 ‘) ~N(t), (22) 
Ms (t) =n(t) —n,(t)-N(t)+n, (t)+n,(t-t, -¢,), (23) 
n, (t)=n,(t-t,)-n,(t-t, -¢,), (24) 
n’(t)=n,n,(t), (25) 
lla =p (8 = a constant), (26) 
dN, 

N(t) =N, (2). (27) 


Equations (19), (21), (24) and (25) have already been explained 
in Case 1. 

Equation (20) states that the number of cells equals the number 
of empty cells plus the number of worker brood plus the number 
of queen brood. 

Equation (22) states that the total number of brood in the nest 
cannot be larger than the number of ‘‘nurses’’ allows. 

Equation (23) states that the total number of eggs which turn out 
to develop into workers and which are laid before time ¢ equals the 
number of occupied cells minus the number of queen brood plus the 
number of workers living at time ¢ plus the number of worker eggs 
laid before time ¢-¢, —¢,. 

Equation (26) states that the ratio of the number of worker eggs 
to the number of queen eggs laid during the same time is constant 
(see Assumption 1). 

Equation (27) states that the number of queen brood equals the 
number of queen eggs laid. This means that we assume that there 
is no egg eating or infant mortality among this brood. This is not 
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true, of course, but we can include these factors by increasing the 
value of the constant B, defined in (26). 

We will now solve equations (19) to (27). We see at once that 
for the interval 7, <¢<T7,+t, the expressions for n(t), 7,(2), 
n,(t), n, (¢) and n,(t) will be the same as those in Case 1. From 
(26) we obtain 


ns(t)=BNQ(t)+y;, (28) 


where y is determined by the known conditions at the time T,. We 
obtain 


y=n,(T,) = Ae*"!. (29) 
From (27) and (28), 
nx(t)=BN(t)+y- (30) 
Equations (23) and (30) give 
(1+ B)N(t) =n(t) 1, (t) +n, (t) +n, (¢-t, -t,)-y 


' & A 
N(t) = No ext _—_— ae iehhin a) nr, eee cy = y 
1+ B Ne No 


for . T, S¢<T, +¢,. 
We have obtained this expression because n(¢), n,(¢), and n, (¢) 
are already known, as mentioned above, and n, (¢-¢, —¢,) will, 
for all the times between 7, and 7, + ¢,, be the number of worker 
eggs laid before some moment in the time interval r,-t,-t, to 
T, —¢,; and this interval occurs earlier than the moment T,. Now 
when we know the expression for N(¢) we can obtain all the other 
functions from equations (19) to (27). We will list the results be- 
low. These expressions are not valid in the entire interval 
T, <t<T, + ¢,, but only in the interval Te GEST ask 


n(t)=n,e**, 


16 


n,(t) = igs +(1- 2)n4| per 
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“ xt 
AU pole 3 
6 


n, (t)=(1- a)n,n,e*', 


na (t) = an, e**, 


n,(t)=n,e** jr a - C=C or) oe rs 
1+8 n,(1+ B) 
A FOE iY 
n,(1+ B) B+i1 
B “f 
ng) = 5 Noe Ja-aa rags ® + 


4 rete) fet ae (30a) 
Bei’ 


No 


N(e) = dd ext E oP a) (1 Za ce) + Ses %) A eat] a 
n 


B+i1 oe Sy 
x 
Bais 
N,@=0, 
N,(¢) = N(é). 


In these equations A, z, and y are defined in (18), (15), and (29). 
The equations are valid even if B =0, that is, when only queen 
eggs are laid. We have not, however, solved the equations (19) to 
(27) for the time interval T, +¢, <¢<T, +t, yet. This we could 
do in the same way as in Case 1 by expressing all the unknown 
functions in terms of n,(¢) and then solving the equation for this 
function, using the known values of n, before time 7, +¢,. How- 
ever, because nobody has ever observed that the time ¢, is larger 
than ¢, in nature (contrary to what is the case among bees) we 
have in the following assumed that ¢,=¢,. This facilitates the 
analysis considerably. If later someone wants to replicate this 
analysis with ¢, #¢,, it will be very simple to do so. It will be 
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more work but the theoretical difficulties will not be greater than 
what they are in this paper. 

The state we have considered here, Case 2, has now given us an 
expression for the number of queen brood existing in the nest. We 
will now proceed to the last case, that is, the time after the first 
reproductive individual has appeared until the time when the colony 
breaks down in the fall. 

Case 3. TT, +t, st <7. The equations for the state of the 
colony will be: 


n,(t) = &n(t), (31) 
n(t) =n, (t)+n,(t)+N(d), (32) 
n,(t)=n,(%) +n, (2); (33) 
n,(t) + N(t) = =e (34) 
n(t) = n(t) —n,(t) —N(t) +n, (t) +n, (t¢-t, -¢,), (35) 
n,(t)=n,(t-t,)-n,(t-t, - ty), (36) 
n’(t) = n,n, (t), (37) 
i -B, (38) 
N,(t)=N,(t -¢,), (39) 
N(t)=No()-N, (4). (40) 


Equations (31) to (38) were explained earlier. 

Equation (39) states that the number of queens at time ¢ equals 
the number of queen-producing eggs laid before time ¢ -¢ 4» Which 
is self-evident by the definition of ¢, and the assumptions we have 
adopted. 

Equation (40) states that the number of queen brood equals the 
number of queen-producing eggs laid minus the number of queens 
which have already appeared. 

Now we will give the solution of equations (31) to (40) as usual 


by expressing everything in terms of nm, and solving for this func- 
tion. From (38) we obtain 


ng (t) = BN, (t) +8. 
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Here 5 can be calculated from the known values of n, (¢) and N, (¢) 


at time he + t,. We find 


O=n,.(T,+t,)-BN,(T, +t, =y. 


That is 
ns(t)=BN,(t)+y. 
From (31), (32), and (34), 
n, (t)=(1-4)n,n(t). 
From (43), (33), and (37) we obtain 
1 


(1 ~ X)ngn(t) =n, (4) -—n'(2). 
6 


Equations (39) and (40) give 
N(t)=N,(¢)-Na(t-¢,)- 
Equations (45) and (42) give 


N(t) = : In, (t)-n,(t-t,]. (6 #0). 


From (46) and (35) we find 


B+il 
B 


Equations (31), (47), and (36) give after a few simplifications 


(1 - a)n(Z) = Pan g(d- : n,(t—t,) —n, (¢-#,)- 


Finally we obtain from (48), (36), and (44) 


1 B+t , erect. a5 CAST ee 
moe Lai np =a) 8 ) 
iS aa(B +X) F 
aioe: he BU Ss ood _— 
masa baht B ng ( 


(n, + t)ng(t~ 0) - Png (t-t,) +n, (2) =#,) = 0. 


n, (¢) - Balt 8) mgt t, — t,) =n(t) —n, (t) + 


(41) 


(42) 


(43) 


(44) 


(45) 


(46) 


(47) 


(48) 


(49) 
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By multiplying all the terms in (49) by the constant 


Ne B (1 > a) 
B+1 
and remembering that we have assumed that 
t,=¢,, (50) 
we can write equation (49) in the following form: 
mi,(t)+(1— a)ngngn, (t) =n, (t-¢,) + 
B(n,+1)-2, 
Te re eo eee oes (51) 


n,B(1- &) 
B+1 


All the expressions on the right side of (51) are known and thus 
we can solve this equation which is an ordinary differential equa- 
tion in n, (¢) of the first order and degree. We can solve for all the 
time interval T, +¢, <¢<T7,+2¢,. Now we can be rather con- 
fident that T lies somewhere in this interval. More than 2¢, units 
of time will probably not pass between the time 7, when the colony 
begins to produce queen eggs and the final breakdown of the colony. 
If, however, we want to continue the solution we can do so step by 
step: (1) For the interval 7, +2¢, <¢<T,+¢,+¢,. (2) For the 
interval T, +¢,+¢, <¢<T,+32,, ete. 

Now we will solve equation (51) for the interval T, +¢, <¢< 
T, +2¢,. By introducing expression (30a) for n, (¢) into (51), we 
obtain 


, B ee 
ns, (t) si non,(1 - X)n, (2) = ER ; n,tBe (tt 4) " 
B(n,+1)-n ae 
Tia Belt ea) Appeen ee eee 
ASE Mac eats | n,B(1— a) (t-t 
soit AF ye No ee = Pag) eT eS x 17%) 
(B + 1)? Me 2y B+1 Ae ) 
where B is defined as 
B=(1-a)(1 # A —x(ty+tg) 
OL eet ant (53) 


oe RM 
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Equation (52) can now be written 


ni(t) +n, (t)=se**+u, 


where 
r=(1-a)n,n,, (54) 
B ~xt B(ng+1)-n 
3= Vo eee ee 1-~a)- 
gee | B+i1 Su) Bs) 
Ang (1 ~ 4) ms] 
n B 
B(n, +1)-1n, 
= —— 1- 
The solution of (52) is now 
n,(t) = ; esto Gere, (57) 
2+r r 


where y is determined by the conditions at time 7, +¢,. We can 
easily see that 


1 SEO AREY [ B ERTIES x 


che ap B +1 
(58) 
8 a ee ee i) 
c+r r 


Now we can obtain the number of queens at the time T, N, (T). 
We obtain from equations (39) and (42) 


N, (2) =Ng(t=t) = 5 ng(t~t,) ~ 5. (59) 


If we now assume that 

ftps T STs ote (60) 
or the equivalent assertion 

tet, Sie PT, 
holds, we find 


“Ty 
n)B gece) _ A 


1 
v,ry~ 3 [myer -y]= 20% Rete (61) 
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This function increases if T, decreases. If we assume that 
T,+2t,<T<T, +3¢,, (62) 


we find that 


- “(T= =stT = 
N, (1) = Abe Slots eri aedl h pea se" oe: (63) 
Bier+r r 


Before we proceed to an investigation of the behavior of N, (¢) 
we will list the other functions in which we are interested. We 
can find all of them from the expression for n,(¢). For T, +t, 
t<T,+2¢t, we have 


B+i 
n(t) = ae In, (¢) -n(t— t,)], 
n,(t) = ta (8 a In, (¢)-—n,(¢-t,)] + 
B+1 y ; 
n,B(1- 4) Ing (42) -ng(¢-¢,)], 


1 
EO eG a 
8 


n, (¢) =n, In, (¢)+N()], 
ng (t) = an(t), 
n,(t)=(1-a)n(t)-N(d), 


8 


u 
n, (¢) = wa a —— es 


e 
+r 


1 
N= 5 -N,O- 4s 


1 
N,(@= gm(t-t))~ 2, 


N,(t)=N(t) +N, (0). 


Now we have to investigate the behavior of the function N , (4) 
for different values of 7, (and 8). We could of course do that 
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analytically but the expressions we obtain in that way will be very 
complicated and cumbersome to work with. We have chosen to 
make a numerical analysis of the form of N,(T) as a function of 7,. 

Determination of the optimal value of T,. We chose the unit of 
time to be one day and we chose the following values of the con- 
stants which occur in expressions (61) and (63): 


¢,=30. This value is given by Richards (1953, p. 49) for 
Vespula. 
t,=40. This value is taken from what we know about bees; 
that ¢, > ¢, for wasps we know with certainty. 
n,=0.2. This seems to be a probable value which has been ar- 
rived at by considering several empirical reports. 
n,=0.7. Richards (19538, p. 68) where he tells that 118 workers 
(Polybia) built 167 cells in two days; perhaps this value 
is higher than what is usual in well-settled colonies. 
= 20. This value is given by Richards (1953, p. 47) for 
Vespula. 
T =150. Approximative length of the summer. 
a =0.2. The value is given for Polybia in Richards (1958, p. 71). 


For 8 we will choose different values. 

From equation (15) we obtain 2 = 0.04606. From equation (18) 
we find A = 21.366. The calculations have been carried out on the 
high-speed electronic computer IBM 650 and in Figure 2 to Figure 
5 have been given the form of N, (T) as a function of 7, for dif- 
ferent values of 6 and in Figure 6 the form of N, (T) as a function 
of B for a certain value of 7,. From Figures 2-5 we can draw sev- 
eral interesting conclusions. 

In Figure 2 (8 = 0.2) we see that there is a very marked optimal 
value of T,, T, o54,=90- If the wasps begin at that time to pro- 
duce queen-developing eggs, with this value of f, they will have 
produced 3353 queens when the colony breaks down. Further we 
can see that if they begin at a time earlier than 7, = 87, then 
N,(T) becomes negative! This seems strange but it means that 
with such a low value of § as we have here, very few worker- 
developing eggs will be laid, and before time 7, + ¢, when the first 
queens appear, most of the old workers will have died and the 
number of workers will have decreased so much that this number 
will not be sufficient to tend the brood. A nest with 6 = 0.2 and 
T, < 87 will end with a catastrophy! We see that if the wasps 
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choose the time 7, just two days before the optimal day they will 
lose about 60% in produced queens but if they are ten days late 
they will only lose about 20%. Such a low value of £ as this will 
probably never occur in nature. 

In Figure 3 8 =1 and there too T, ,,,, = 90. Here, however, 
N,(T) will never become negative. The wasps will lose 10% if 
they begin 2 days before T, ,., or 6 days after. Even this curve 
seems to be too steep to conform with reality. 


N, (7) 


3000 
2000 


1000 


-!000 


- 2000 


-3000 


FIGURE 2 


In Figure 4, however, where B = 5, we have a curve of quite 
another character than in the two earlier diagrams. Here N, (T) is 
steadily decreasing in all the range T, = 60 to 7, = 120. However, 
this decrease is very slow up to T, = 90 where suddenly the value 
of N,(T) begins to fall very fast. In practice the wasps can choose 
any time before 7, = 90 with neither gain nor loss. However, as 
it is no use to have a lot of queens hanging around longer than 
necessary, even here T, ,,,= 90. For this value of T,, N,(T) = 
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671 which is a good, and probably not an unusual, value for a 
successful nest. 

Figures 2-4 exemplify the general fact that the first derivative 
of N,(T) with respect to 7, will always have a discontinuity at 
T, = T -2t, =90 depending upon the fact that N, (T) is defined 
by different sets of equations before and after this value of 7,. In 
Figure 9 this discontinuity still exists at T, = 90, but it is not so 
prominent because the maximum occurs earlier. 

Here perhaps is the time for an important observation. Because 
of obvious biological reasons it is in the interest of the queens 
that they shall be hatched as late as possible. They shall only 
have time enough to be fertilized and then it should be time for 
them to enter their winter sleep. Now for the 8-values which prob- 
ably occur in nature we have seen that the value of N, (T) asa 
function of 7, decreases very slowly before T,=T, ..,.- For B 
between 10 and 30 for example, the wasps would certainly gain by 
using a value of T, less than 90. If they begin too early, however, 
many queens will have to wait a long time before entering their 
winter sleep. There will evidently be an optimal value of T, which 
will perhaps differ slightly from the value obtained according to 
the theory laid down here. This value will then be less than the 
value given by our analysis. 
dN, (T) 


This remark is valid only when ( <0, that is, for 
T,=90 


i 
B > 4 or 5. We do not know if values of 8 larger than 5 occur in 
nature but it is possible and perhaps in unsuccessful nests probable. 

In Figure 5 6 = 30 and here the value of N, (T) is decreasing 
rapidly for all the range of possible T,-values. 

Because T, = 90 seems to be an optimal value for the B-values 
which are most likely to occur in nature we have chosen this value 
of T, and plotted V,(T) as a function of 8 in Figure 6. We see that 
the resulting curve is a hyberbola with the codrdinate~axis parallel 
to the asymptotes. The range of B-values that occur in nature will 
certainly be included in the range shown in Figure 6. The limiting 
values of 8B shown in Figure 6 correspond to 200 and 1000 as 
limiting values of N, (7). 

From Figure 6 we can see that the constant B in a certain sense 
measures the success of the colony of wasps. From this figure the 
immediate conclusion is that the nest with the least value of 6 is 
the most successful. Now, however, the wasps cannot change the 
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value of 8 as they want. We do not know if it is known what de- 
termines the value of 8 but the rate with which the queen(s) in the 
nest can lay eggs is of course an important factor. Besides, we 
see that for small values of 6 the peak of the curve N, (T) versus 
T, is very narrow. Therefore, if the wasps do not choose the value 
of T, within one or two days from 7, ,., they will not gain any- 
thing but perhaps all the colony will fail and end with no queens 
at all. As the summer certainly does not last for exactly 150 days 
it would be very dangerous to try to have B too small. We can, 
however, state that the colony which is likely to be most success- 
ful is that which tries to make £8 as small as possible within the 
interval shown in Figure 6. 

We will now choose 7, = 90 and 8 =5 because these values 
seem to be highly probable in nature and compute the values of 
some other interesting functions. 


FIGURE 5 
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We obtain for example (T = 150), 


n(T) = 15768, n,(T) = 3474, 

nx (T) = 15218, N, (fT) = 671, 

n,(T)/n(T) = 0.22, n,(T)/n, (T) = 0.38. 
n,(T) 


Sa Oe 
n,(T)+N(T) 
Richards (1953, p. 70) gives the following observed values for 
Polybia: 
1st nest: n = 7758, 
nm, = 1400, 
andso n,/n= 0.18. 
2nd nest: n = 21600, 
n, = 6994 (plus 93 egg-laying queens), 
andso n,/n = 0.82. 
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Our values of n and n, are between the values in these two cases 
and so is the value of n,/n. As seen in Figure 8, according to the 
theory presented here, the value of n,/n increases for increasing 
time (that is, for increasing n). The values given by Richards are 
thus in good agreement with our theoretical values. 

In Figure 7 1, (¢) and N,(¢) are plotted as functions of time ¢. 
The time is measured with one day as unit and on the vertical 
axis the number of insects is given. For example, if we had plotted 
n,(t) against ¢ for T, = 86 and B = 0.2, that is, for values of 6 and 
T, for which N, (T) <0 in Figure 2, we would have arrived at the 
same curve as in Figure 7 up to ¢ = 120, but between ¢ = 120 and 
t = 150 it would have looked like the dashed curve in Figure 7. 

Finally, for 8 = 5, T, = 90, we have plotted against time a few 
other quantities of interest (Figure 8). We can see here that after 
the time ¢ = 120 the number of workers per cell as well as the num- 
ber of workers per egg, larva, and pupa and the number of ‘‘builders’’ 
per ‘‘nurse’’ steadily increases. Before ¢ = 120 all these quanti- 
ties are constant. That is not a weak point of this theory but de- 
pends on the explicitly stated assumptions in the beginning of 


3000 


2000 


FIGURE 7 


142 BIRGER LOVGREN 


° 100 120 140 160 t 


FIGURE 8 


this paper. If we change the Assumption 5, that is, we assume 
some functional form of n, other than a constant, then none of the 
quantities in Figure 8 will remain a constant for ¢<120. This 
would give us a very handy method of testing Richards’ (1953, p. 
71) hypothesis that the sudden start of queen production depends 
on the fact that the proportion between the ‘‘nurses’’ and the num- 
ber of brood increases (and, as a consequence, the brood would 
obtain better food, or some similar advantage). The quantity n, 
certainly is very difficult to measure. The theory presented here 
allows us, however, to obtain the functional form of n, in an in- 
direct way. We just measure one of the quantities n(¢), n, (¢), or 
n,(t) every day (which also is a difficult task, indeed) and as- 
suming values for 4, n,, no, ¢,, and ¢, calculate the value of n,. 
Or if the value of 4 is dubious we measure two quantities, for ex- 
ample, n(¢) and n,(¢) and calculate n, and «. We can do the same 
even for three constants. We will then see if n, increases and if 
the hypothesis is correct. 

To check if the wasps choose the optimal time to begin the 
queen-production should be relatively simple. We have only to 
measure the constants of the nest and the time they begin the 
production and see if the values are in agreement with those ob- 
tained theoretically here. 

If someone for some reason wants to change to some other func- 
tional form one of the quantities assumed constant here, this is 
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relatively simple and the procedure used in this paper does not 
need to be changed. However, even now the actual computations 
are very laborious and time-consuming, and it would have been 
very cumbersome to make them by hand. We are sure that it would 
suffice to replace the constants here by piecewise constant func- 
tions. Then we can use the results obtained here directly just by 
changing the values of the constants as we want now and then. 

As mentioned above, the computations required are simple in 
principle but very laborious. All the numerical results in this 
paper have been obtained by using the IBM 650 computer at Cornell 
University. 


PART 2 


We will now continue our mathematical description of social 

wasps by treating the case when they build special cells for the 
queens. We make first: 
Assumption 6. The queen-developing egg will be laid at the same 
moment the workers begin to build the queen cell. The error we 
make will be negligible compared with the time it takes for a 
queen to go through all the egg, larva, and pupa stages. 

We will use the same symbols as in the previous part plus a few 
new ones: 


n(t) is the number of worker cells at time ¢. 
N, (t) is the number of queen cells at time ¢. 
N. is the number of queen cells a worker can build per unit of 
time. N, is a constant. 
a(t) is the number of ‘‘builders’’ building queen cells, living at 
time 7¢. 


6 


In this part we will not give any detailed derivations but only 
state the equations and the results. The derivations are almost 
identical with those in Part 1. 

Case 1. 0<t<T,. We will obtain exactly the same values of 
the functions of the colony as in Case 1 in the previous part. In 
addition we will have: N, (¢) = 0, z(¢) = 0. 

Case 2. T,<t<T,+t,. The equations for the state of the 
colony are: 


n,(t) = an(t), (64) 
n(t)=n,(t)+n,(2), ; (65) 
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N,(t)=N,(¢), (66) 
N,(t)=N(%), (67) 
n,(t)=n,(¢) +7, (¢), (68) 
pp Oren (0) mae os (69) 
ns 

nx(t)=n(t)—n,(t) +n, (+n (t-%, -t,); (70) 
n,(t) =n, (¢-¢,)-n,(¢-—¢, -—¢,), (71) 
oti =a (a = constant) , (72) 
dn 

aN , (t 

30 Neel, (73) 
fe =ngln, (t) — x(?)]. (74) 


Equations (64), (65), (68), (69), (70), and (71) were already 
explained. 

Equation (66) states that the number of queen cells will equal 
the number of queen eggs laid up to the time ¢. This depends on 
Assumption 6 above. 

Equation (67) states that the number of queen cells equals the 
number of queen brood (Assumption 6). 

Equation (72) is the mathematical formulation of Assumption 2 
in the Introduction. 

Equation (73) states that there cannot be built more queen cells 
than there are workers available for the job. 

Equation (74) states the same thing for the worker cells. 

In the solution we proceed as usual. The equation for n, (¢) 
becomes: 


ngN,(1- a+ 4) 


P No P 
ne (t) +n, (0) ant (b=8g 


N,+7,@ 
n,N,(l1-AK+n, -n, Xk +n, 
nities ee A) ‘i (75) 
8 6 
nmgN,.(1-4) n,n, dN,(1-0 
y(t time) oten’ owes ie eae 


N,t+nga No+n,a@ 
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where 


Ty 


b= -an,e* (76) 


After solving equation (75) in the usual way we can obtain all the 
other functions in which we are interested. 


Case 3. T,+t,<t<T. The equations for this last stage in 
the development of the colony are: 


n, (t) = an(t), (77) 
n(t)=n,(t)+n,(¢), (78) 
n,()=n, (+n, (4), (79) 
n,(t)+N(t)= Salt) (80) 
Ns 

n,(t)=n(t)-—n,(t)+n, (+n, (¢-¢, - 2), (81) 
n,(é)=n,(¢-¢%,)-n,(¢-¢, -¢,), (82) 
fae =a (83) 
dn ; 
aN ,(t 
sei =N,2(2), (84) 
dn(t) 

ae n(n, (t) - x(¢)], (85) 
N,(t)=Ng(t-t,), (86) 
N(t)=N,(Q-N, (2), (87) 
N, (4) =N, (2) +N, (2) (88) 


By expressing all quantities in n,(¢) and so finding the equation 
for this function we find that it is identical with equation (75). 
From (86), (87), and (88) we find 


N,()-N,(¢-t,) =N,(t-#,)- (89) 


If a@ = ~, which it is in most cases in nature, equation (75) is 
not valid. We will indicate the solution even in this case. We as- 
sume as in the earlier part that ¢, =7¢,. For Case 2, that is, when 
T, <¢<T,+,, we obtain the results: 
x 
n(t)=nge }, (90) 
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n,(t) = an,e*'}, (91) 
n.()a(1- ange” *, (92) 
n,(t) = de* 1 4 (1-a)n,e*'!, (93) 
n(=[2 + - aa 1, e**, (94) 
6 
Ng (t)= —— e*#- — + pe", fod. 
e+p p 
where 
p =nNg; (96) 
2 97) 
o=N, peru Sete) UT Nos ( 
6 
r=N,n,(1-a)n,e" 1, (98) 
es ent : en . (99) 
p @+p 


In Case 8 the expressions for the functions of the nest are rather 
complicated, and we will only give the expression for N,(¢). That 
is, for T, + ¢, <¢ < T we have 


) 
Na(fim gam ett ner, Gon) 
where 
6=AN,e “\1(e°*'! — e7*'8) (101) 
b= Ngno(t-a)(1~n,)e*™, Oak 
gh a i pee ° Poabd Wht cae poe POY cet heb) be (103) 
w+ ‘ 
p P 


From (89) we can finally find the expression for N, (t). We will 
now make a numerical analysis of N, (¢) using the same values of 
the constants involved as in Part 1. For the new constant N g we 
assume the value 0.5 quite arbitrary, just because the queen cells 


are larger than the worker cells. We will make the analysis for 
a@=oo, Thus N,=0.5, a= 0, 
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In Figure 9 N,(T) is plotted as a function of T,- We see that 
T, opt, = 86. This means that wasps which build special cells for 
the queens should begin the queen production earlier than wasps 
which use the same cells for all individuals. 

The number of cells, wasps, etc., can now be determined for 
every time. Due to lack of space we cannot give any further de- 


N, (7) 
3000: 


60 70 80 90 100 10 120 
FIGURE 9 


tails but as the procedure is the same as in Part 1 above anybody 
can carry out the necessary numerical calculations. The same re- 
marks and conclusions as at the end of Part 1 are also valid here. 
Several hypotheses about what time the wasps choose to begin the 
production of queens can be checked by the theory given here. If 
someone wants to change the constants in this paper this is per- 
fectly possible. There is particularly one constant, 4, which it is 
perhaps necessary to change. After the time 7, it seems probable 
that there are not so many worker eggs laid as before that time, 
This is equivalent to an increment of % because % is the propor- 
tion of the worker cells which are empty. If we assume some 
simple form of the time dependence of «, the procedure above will 
be slightly changed but hardly more complicated. 
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In all this paper we have only considered queens and workers. 
However, males will also be produced. It has been tacitly implied, 
however, that we have classified them together with the workers. 
Now of course they will not work but the error we make in this way 
is negligible because the number of males is very small compared 
with the number of workers. 

In this paper the main importance has been attached to obtaining 
and explaining the mathematical theory developed. The numerical 
calculations have only been performed in order to illustrate, il- 
luminate, and in a certain sense justify the mathematical theory. 
Therefore the values of the constants have been chosen not for one 
certain species but as a kind of average values for all species of 
wasps. Ihe values have been taken from Richards’ 1953 because 
this book was easily accessible to the author during the writing of 
this paper while Richards’ 1951 was not. It is the intention of the 
author to use the detailed accounts of different wasp societies in 
Richards’ 1951 (and in other sources) in order to make a careful 
comparison between the observed values given there and the 
theoretical values derived in this paper. This will be done ina 
later paper. 
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Using records of the response of the extraocular muscle system to a 
step input, the frequency response characteristics of the system have 
been calculated. They show a progressive reduction of gain for fre- 
quencies above 30-40 cps. The high frequency tremor of the eyes ob- 
served during ‘‘steady’’ fixation would thus be expected to be damped 
down and this is in fact borne out by frequency analyses of tremor records. 


An eye movement response to the sudden displacement of a 
stimulus (light spot) from one position in the visual field to another 
is illustrated in Figure 1. There is a reaction time of the order of 
0.12-0.18 seconds. Enough evidence (Monnier, 1952) is available 
to permit the conclusion that most of this delay occurs at the 
sensory and central levels and that it may be ignored in analyses 
involving only the motor system of the eyes. In attempts to de- 
scribe the feedback mechanism between the retinal images and the 
motor system, this delay, which is somewhat variable, cannot, of 
course, be neglected. Once initiated, however, an oculomotor re- 
sponse has predictable characteristics; only its time of onset 
needs a probabilistic treatment. 


s 


aor ene 


FIGURE 1. Movement response of the two eyes following sudden dis- 
placement of a light spot in the visual field through 20°. ‘Time line in- 
terrupted every 1/100 sec. Stimulus change at S. Positions of eyes at 
beginning and end of record corresponds to 20° movement. 
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The evidence (Westheimer, 1954a; Breinin, 1955) overwhelmingly 
favors the formulation that the response illustrated in Figure 1 re- 
sults from a simultaneous change-over of innervation to all the 
extraocular muscles, i.e., from a step input to the extraocular 
muscle system. 

It has been postulated (Westheimer, 1954a) that as a first ap- 
proximation the system is linear and satisfies the equation 


A,6+A4,6+4,9=f(0), (1) 


where @ is the output (i.e., the eye position in the orbit), the A’s 
are constants, and f(t) is the forcing function (a function of time). 
For simplicity this discussion is confined to a unidimensional 
situation. 

Equation (1) may be rewritten, in conformity with the notation 
adopted by many authors (Trimmer, 1950), 


6+260,0+0720=? f(t), (2) 


1 
where w? =A,/A, and €=A,/2(4,A,)?. The parameters ¢ and 
w, are descriptive of the system. 

If Figure 1 is accepted as a response of the system to a step 
input it is possible to evaluate the parameters. Representative 
values for the human extraocular system are ¢ = 0.7 and w, = 240. 

Another way of studying systems characteristics is to investi- 
gate the system’s response to sinusoidal forcing at various fre- 
quencies. This, for example, has been employed recently in a 
servoanalytic study of the pupillary reflex (Stark and Sherman, 
1957). Such an analysis is not readily feasible in the case of the 
extraocular muscle system. Eye movement responses to visual 
stimuli usually have time characteristics leading to the conclusion 
that they result from inputs to the extraocular muscle system which 
are either step functions or step velocity functions, producing 
saccadic and smooth following movements respectively (Westheimer,: 
1954b), 

The parameters obtained from the response to a step input al- 
low us, however, to predict the amplitude and phase relationships 
between the eye-position response and the sinusoidal input (Trim- 
mer, loc, cit.). They are shown graphically in Figure 2. For fre- 
quencies of several cycles per second (cps) no significant dif- 
ferences are expected between input and output, but above 40 cps 
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a given input amplitude gives rise to a smaller response amplitude 
and a phase lag is seen; these effects increase with increasing 
frequency. 

The fastest known gross sinusoidal eye movements occur in 
voluntary nystagmus where oscillations of 6° amplitude and fre- 
quency of the order of 10 cps may be observed (Westheimer, 1954c). 
This is still within the region in which no significant differences 
between input and response are expected. It follows that when 
considering gross eye movements, as for example in studying the 
feedback characteristics between the retina and the oculomotor 
system, one may regard the response characteristics of the latter 
as practically ideal. 

Recently considerable attention has been given to the small eye 
movements that occur even during ‘‘steady’’ fixation. Several 
clearly distinguishable components of this micronystagmus have 
been described; none exceed a few minutes’ of arc in extent. The 


AMPLITUDE 


PHASE 


AMPLITUDE RATIO (OUTPUT/INPUT) 


-90 


PHASE (DEGREES) 
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FicureE 2. Theoretical response characteristics of extraocular muscle 
system to sinusoidal forcing, based on parameters deduced from the re- 
sponse to a step stimulus illustrated in Figure 1. Ordinates: ?—Ampli- 
tude ratio (response/input), L—phase lag (degrees). Abscissae: fre- 
quency of sinusoidal input (cycles per second), 
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one component of interest in connection with this analysis is the 
fine sinusoidal tremor (Riggs and Ratliff, 1951; Barlow, 1952; 
Ditchburn and Ginsborg, 1953) that has an amplitude of only a 
fraction of a minute of arc and a frequency range from 30 cps to as 
high as 150-200 cps. Figure 2 shows that at the higher frequencies 
the input strengths would have to be relatively enormous or the re- 
sponse amplitudes quite small. 

The results of Riggs and co-workers (1954) indicate a reduction 
of tremor amplitude when the frequency is high. Fender (1956) 
carried out a harmonic analysis of his tremor records and found a 
‘spectrum extending from 30 to 150 cps with a maximum at 35 cps.”’ 
Hebbard (personal communication) using a method of inspection 
measured the amplitudes of several thousand individual cycles of 
tremor and found that the highest maximum amplitudes occur when 
the tremor frequency is 30 to 50 cps and that the maximum ampli- 
tude at each frequency decreases as the frequency increases, 
tremor cycles of frequencies of 140 cps and up having amplitudes 
near the resolving power of his apparatus (2-4 secs of arc). 

A classical dictum in the design of servomechanisms is that a 
compromise must be struck between faithfully following the input 
signal and ignoring the noise. In the case of the extraocular muscle 
system the input consists of impulses from the supranuclear cen- 
ters in the central nervous system which control eye position. 
Discrete nonsynchronous firing of motoneurons is seen in electro- 
myograms (Breinin, Joc. cit.) to produce high frequency activity 
with components extending into the range of several hundred cps. 
Peaks separated by intervals of the order of 10 msec or less seem 
to be particularly prominent. The consequent nonsynchronous con- 
traction of the muscle fibers must inevitably result in fluctuations 
in muscle tension and it is here that we may look for the origin of 
the tremor. 

The fact that the electromyograms have quite strong components 
in the frequency region in which the tremor amplitude progressively 
decreases constitutes a verification of the sinusoidal response 
characteristics deduced from the step input response and shown in 
Figure 2. In damping down the high frequency fluctuations the 
oculomotor system treats them just as noise is treated by a well- 
designed servomechanism. 
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The spontaneous activity which is common to several natural afferent 
systems (thermal receptors, semicircular canal receptors) may be sup- 
posed to reduce the sensitivity of detection of stimuli by the organism. 
The simple device of reciprocal inhibition (conjoined negation) can be 
shown to suppress spontaneous activity. To illustrate this principle, a 
McCulloch-Pitts network is constructed for thermal sensitivity and its 
mathematical properties are examined. 


It has been shown that each hair cell on the rabbit ear can acti-. 
vate on the average about five ganglion axons, but each axon can 
be activated by about 80 hair cells (Weddell, Taylor, and Williams, 
1955). There are about 100,000 hair cells to about 5,000 ganglion 
fibers. It has been shown (Williams, 1958) that if the cortical 
interconnections are similar to the peripheral interconnections then 
this multiple innervation limits the effect of chance axonal failure. 
The central pattern which is produced by stimulation of a single 
hair is only slightly changed when an axon fails. However, this 
pattern may be affected by an accidental firing of an axon. It is the 
purpose of this note to discuss this effect. 

The effect of chance firing of axons in the hair system can be 
considered with respect to 1) localization and 2) cause. 

1). A single accidental firing of an axon would produce not more 
than about 80 individual central activations. These diffuse activa- 
tions are not adequate for accurate.localization. Thus no single 
hair can be localized by a single chance axonal firing. In this 
case, error is not amplified with respect to localization because 
the (false) information generated by a single accidental firing is 
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dispersed in a pattern which is meaningless by itself. This con- 
ception is to be contrasted with a model system in which each hair 
is linked by a single axon to a single central cell without exterior 
connections: chance firing and chance non-firing are reproduced 
absolutely. 

2) On the other hand, information with respect to cause is am- 
plified because any activation at all in such a system could only 
have been caused by mechanical deformation of the hairs. How- 
ever, it is questionable how much meaning is conveyed by a pat 
tern which is indeterminate with respect to localization and fre- 
quency, even though cause be predetermined. 


FIGURE 1 


In the case of spontaneous activity where frequency may be of 
an appreciable magnitude, the problem of error control with respect 
to cause may become more serious. R. Granit (1955) has dis- 
cussed spontaneous activity in several systems, and although 
there are many features which these systems have in common, this 
note will consider only the case of thermal receptors. 

Y. Zotterman (1953) distinguishes between two classes of spon- 
taneously active thermal receptors in the tongue of the cat: (1) The 
activity of the first category is increased by cooling the surface 
of the tongue and is inhibited by warming; these are called ‘‘cold 
receptors.’’ (2) The activity of the second category is increased 
by warming the surface of the tongue and is inhibited by cooling; 
these are called ‘‘warm receptors.’’ The construction of a simple 
McCulloch-Pitts net for these two classes of fibers using only 
reciprocal inhibition (conjoined negation) provides a simple effec- 
tive method for the suppression of spontaneous firing, irrespective 
of whether the firing is truly random or is caused by non-iso- 
thermal conditions across the terminals (Lele e¢ al, 1954; Tyrrell 
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et al, 1954. See Figure 1). It is reasonable to suppose that the 
two systems could be joined at the several levels of integration 
(dorsal columns, thalamus) by such a mechanism. A further con: 
sequence of this net is that an increase of activity in either one or 
the other channel immediately clears the opposite channel of all 
activity for a period of time which we shall call 5. 

If the mechanism of Figure 1 is symmetric and has identical 
elements with no variability in delay time then the higher levels 
of cross inhibitory elements would introduce additional effects. 
The frequency of firing in each chain would be reduced at the first 
level only, the reduction being the same for both. If, however, 
there is an element of randomness introduced along the chain the 
result is different. If there is statistical independence at each 
step we can calculate the result. Statistical independence at each 
step can be approached most easily if we suppose that the mech- 
anism of Figure 1 represents only the most probable interconnec- 
tions, the fibers actually going different distances along the chain, 
these distances being randomly distributed. In the figure the aver- 
age number of inhibitory per excitatory fiber is one. If it is some 
larger value, R, we need only replace 6 by a value very nearly 
equal to fF § in the subsequent equations. 

From the network shown in Figure 1, we next calculate the fre- 
quency v,, of the output of a, as a function of the frequencies v, 
and v, of the afferent neurons a and b, All thresholds are assumed 
to be equal to two. Assuming statistical independence of the fir- 
ings of a and 3, for the reasons given above, we find (Landahl, 
McCulloch and Pitts, 1943). 


Ve, = (1 -»,), (1) 
Va citiy = Yaill — 85,4)» (2) 


where 8 is the duration of the inhibitory effect. Similar expres- 
sions can be found for the frequencies of the 6 chain by inter- 
changing the subscripts a and 6 in the above expressions. By 
successive substitution we find the following expressions: 


v =v, —25v, 4 + 82 Vy (Yq +4) — 89 (Ma 4B)» (3) 


a2 


2 fe 
Vag = Ug ~ 89UQY, + 357, v, (vy + v5) 


(4) 


2 
devi v,(v2 + Ty vy tp) teres 
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2 
Vag =Uq ~ £5vgv, + 8S viv,(v, +) - 


(5) 


25° viv,(2ue + llviv, + Que) +... 


Because of the symmetry in expressions (3) to (5), v,; - v, must 
equal v,, — v, so that 


A cna, hee el (6) 


Since v,, and v,, both decrease with 7, it can be seen that, if 
v, > v,, the average frequency of a, will approach vy, - v, while 
that of 6, will approach zero if the chain is sufficiently long. 
Since the above expressions are cumbersome it is useful to ob- 
tain approximate expressions for the frequencies. Because of (2) 
and (6), we may write 


Va, = Var Ye Fe? (7) 
A Vas = Macinty ~ Ya = ~ OM aiMog = — OMG; (Vi — Ma + %)- (8) 
Since Ai = (¢ + 1)-i=1, and since Av,,/Ai =dv, ,/di for small 
values of 5, we have approximately 
avy; 
di 


Since v,, =v, for ¢=0, the output of a, is found from (9) to be 
given by 


> bv; (Ya; = Ngot vy)" (9) 


. Va (v4 is vs) 0 

eon ye ye taker WB? (1 ) 

The output of 6, can then be found from expressions (7) and (10). 
Note that when n5(v, — v,) is small, we have 


= 1 53,2 
Mag ue —ndv vz, +.— 06 n Var, (vy +R) — eee (11) 


For n=, vj, =v, —v,+ The output of chain a is then approxi- 
mately the difference between the input frequencies while that of 
chain } is approximately zero. 

To find the length of the chain necessary for vi, to approach 
the asymptotic value v, - v, within a certain fraction E of that 
value, we put v,, —(v, — v4) equal to E(v, - v,) and solve for n. 
This particular value, n*, of n is found to be given by 
1 (1-E)v, 


= —————_ lo 
50,-%) 8 Ey, 


n* 


(12) 
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From this expression it follows that the smaller the difference be- 
tween the input frequencies, the longer must be the chain in order 
that only one side should have an appreciable output frequency, 
On the other hand, the longer the inhibitory period §, the shorter 
the chain required, But if 5 is fairly large, the approximation of 
the difference equation (8) by the differential equation (9) is no 
longer satisfactory, If 6 is not small enough the chain length re- 
quired by (10) is actually longer than required by the solution of 
(8). This can be seen from the fact that the second term in the 
expansion of (10), given in (11), contains the numerical coeffi- 
cient n?/2 which is larger than the corresponding numerical co- 
efficient in expressions (3) to (5), But since the difference equa- 
tion (8) also breaks down for too large values of 5, we shall not 
pursue this problem further, 

In afferents from the semicircular canal system spontaneous ac- 
tivity occurs bilaterally. Movement in a particular direction in- 
creases activity in one side and inhibits the activity in the op- 
posite side. If these systems were joined by a similar mechanism 
to the one described above, the spontaneous activity would not be 
expected to reach higher centers. The model which we have de- 
scribed may therefore be common to many systems, It is also pos- 
sible that a similar device might be constructed for the suppres- 
sion of ‘‘noise’’ in communication systems. 


This work was aided in part by a grant from the Dr. Wallace C. 
and Clara A, Abbott Memorial Fund of The University of Chicago. 
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On the basis of a simple model and by introducing rather natural as- 
sumptions, it is shown that it may be possible to estimate relative blood 
flow rates if blood pressure values are obtained simultaneously with 
electroplethysmographic readings and if several sets of measurements 
can be made under suitable experimental conditions. 


It is the purpose of this note to consider the possibility of esti- 
mating cerebral blood flow from electroplethysmographic measure- 
ments (cf. Nyboer, 1944, 1950; Bonjer, 1952). We shall attempt to 
do so in terms of the simplest possible model of the cerebral blood 
vessels. Only that part of the vessel which is very distensible 
and which shows marked volume changes is of interest. There is 
evidence to indicate that attention should be focused on the ar- 
terioles (cf. Hoerr, 1944, p. 1383). We shall suppose that the blood 
enters from the small arteries where the pressure is P(t) (Figure 
1). The mean pressure in the distensible elements is denoted by 
P,(t) while that at the venous side is denoted by P,,(2), this latter 


| a a a aS ah a 


P(t) —> P. (t) —> Rit) 
SEES oa Be Bae re aes se ee ee ec ke Sk AEE: 
FIGURE 1 


being approximately pulse free and ordinarily fairly small com- 


pared with P,. 

Denote by p the sum of the resistance of the small arteries and 
of one-half of the resistance of the distensible element. Similarly 
denote by o the sum of the resistance of the capillaries and ven- 
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ules, and of one-half of the resistance of the distensible element. 
Denote by A the volume of the distensible element, This volume 
changes with time when ‘there is a difference between the inflow 
and outflow. The inflow is equal to the pressure difference P, — P, 
divided by the resistance p while the outflow is given by the pres- 
sure difference P, - P, divided by the resistance o, Hence we ob- 
tain the following differential equation: 


A t i v v 

——= + - += ++ “+ - —P., 1 

dt p g a p po () 
Next we shall suppose that the distensibility is such that the 

volume A and the internal pressure P, are related according to the 

expression 


A=A,P,+K, (2) 
where A, and K are constants. Substituting (2) into (1) we find 
FAS: see 
—— a a+ -—3 = BA+ BR, (3) 
dt a p 
where 
ptoa 
= . 4 
aires (4) 


Both p and o are functions of A, since the resistance of the ele- 
ment depends on its volume. However, we shall neglect this ef- 
fect and consider constant average values of p andc. 

We assume that the electric conductivity C across the brain is 
increased in proportion to the increase in the amount of fluid 
in the small blood vessels distributed throughout the tissue, 
We shall also assume that the number of active capillaries is 
constant (cf. Hoerr, 1944, p, 182), Denoting by the subscript 
m a normal mean value for a given individual, and by C, a con- 
stant, we may write approximately 


C=C +Aa(A-A )=C,+aA, (5) 
Putting 
o=(C-C,), A=O/Bp, (6) 


and denoting by P, the mean value of P, , we obtain from (5) and 


(3) 
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t 
o(t) = A pe-® eft’ (P, — P,) dt. (7) 
As an example of how the parameters can be evaluated, consider 
that P, is a periodic, sinusoidal pressure wave with an angular 
frequency 27/T. 
Actually T should be about 4 second when the heart rate is about 
1/sec. If P| is the amplitude of the pressure then (7) shows that 


ABP, 


- 2? 7? 
BY + 73 


c(t) = sin — (t—r), (8) 


so that the conductivity lags by a time 7. If r<<T, it can be 
shown that 


B =1/r. (9) 
Suppose that C ~ C, and that the value of C, is not very different 
for various individuals. Then the amplitude c, of the conductivity 
c, is given, according to (8), by 


ABP a P 
9 eH Ri geee h~= bFgs (10) 


, i iy 


The quantity c, is the measured amplitude of the conductivity 
change, P_ is the observed amplitude of the pressure, 7 is the 
period, and B is known from (9). Hence %/p can be estimated from 
the amplitude of the conductivity: 


(11) 


We next consider the case (experimental situation I) where it is 
possible to block the arterial blood flow for a short time. In this 
case, since the inflow is blocked, there is only an outflow term. 
Instead of expression (1),we now have 


Ei~ P. 
de GA,P, + ak — aA, | c (13) 
dt oA, oA, 


because there is no inflow or outflow via the arterioles. If the 
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arterial blood supply can be stopped long enough to estimate the 
time constant then oA, can be estimated from (18) and hence 
knowing £, pA, can also be calculated from expression (4). From 
the product of &/p and pA, we obtain aA,. If the arterial supply 
can be stopped long enough to enable also the estimation of the 
asymptotic level, then we have 4A, P,+ 4K — aA, and hence we 
may estimate the quantity a(K - A_). 

Another way to obtain additional information on %/o (Experi- 
mental Situation II) is to consider an abrupt change in internal 
pressure, as may be produced when the pressure in the lungs is 
suddenly increased. The mean value of the conductivity would be 
expected to change. Denote by primes the values of parameters 
and variables when the pressure is increased by an amount 6P. 
Thus P’=P,+65P. The arterial pressure may be altered but we 
suppose that it is measured when necessary. If the parameters B 
and &/p change we suppose that they are measured under the new 
situation. If C is the mean value of the conductivity then the 
average change will be 


eee (Sr-a(Jroer-e(2)e an 


If P’ is substantially larger than Bs we ignore the change in (d/o) 
which may arise from the change in P,. Then all quantities in (14) 
except 4/o can be measured, and thus 4/o can be estimated from 
C’ —C using (14). The ratio of «/o to &/p is p/o so that the rela- 
tive importance of p and o can be estimated from this situation. 

In addition to the above procedures, it may be possible to block 
the outflow momentarily. This case can be treated in a manner 
similar to that of the first experimental situation. However, since 
P, Varies markedly with time, the interpretation is not as simple 
in this case, 

Since we are here investigating only the possibility of measur- 
ing cerebral blood flow on a relative basis, therefore we consider 
flow F* under given conditions in relation to the mean flow F 
under ‘‘standard’’ conditions. From (1) we have 

ES Col Tt SEP Coe CE Ce NS 


_—_ = =—-_—- ) ———————_- = 


F P,-P, ot +p* Py =P .(cA,)* +(pA,)* e 


v 


We have assumed that the total number of capillaries is not 
changed. We also have assumed that P,, and hence FP, , can be 
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measured. Since P,,is relatively small, we either ignore it or use 
an estimated value for it, ignoring differences between P* and P . 
We then consider only the term (p* + o*)/(p +c). If ‘ is ove 
siderably larger than o, as determined from considerations dis- 
cussed in the last paragraph (Experimental Situation II), then 
(p* + o*)/(p + oc) = (X/p)/(4/p)*. If data are available for both 
situations, then we can calculate (oA\,); (eA,)*;-(A,), (p A;)* 
and A,/A* (from 4A ,/aA*). Thus we find the relative flow rates. 
It can be seen that it is plausible to expect to be able to estimate 
relative flow rates from electroplethysmographic and arterial pres- 
sure measurements. Even though it is not convenient or possible 
to make all the necessary measurements, it may be possible, by 
correlating the results with other methods, to make approximations 
which will yield rapidly sufficiently accurate results. The pos- 
sibility of making absolute estimates of blood flow is less en- 
couraging, though if it is possible to estimate « independently 
this may be feasible. 

We have considered only the case of cerebral blood flow where 
the volume is indirectly estimated. The methods suggested here 
may be applied wherever volume measurements can be made di- 
rectly (Perry, Hawthorne, Kadetz, and Marbarger; 1957). In terms 
of the above model this means that we let C be volume so that 
C, = 0 and & = 1 throughout. In this case the results will depend 
on how satisfactory are the approximations. If they are satisfac- 
tory and if a sufficient number of different measurements can be 
made then it becomes possible to estimate the average rate of 
flow. 


The author is indebted to Professor T. Postelli for raising this 
problem and to Dr. Robert Macey for several valuable suggestions. 

This work was aided by a grant from the Dr. Wallace C. and 
Clara A. Abbott Memorial Fund of The University of Chicago. 
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As suggested in previous publications, freedom may be defined quanti- 
tatively as a restriction upon the choice of a number of activities. If the 
choice is determined by maximizing the satisfaction function, it is sug- 
gested that freedom may be defined in terms of the satisfaction function. 
If an individual is isolated and no physical restrictions limit his choice 
of activities, he is free to choose any activity in an amount which maxi- 
mizes his satisfaction. This isolated state may be considered therefore 
as that of maximum freedom. If the individual interacts with another, he 
will choose different amounts of his object of satisfaction depending on 
whether he behaves egoistically or altruistically. But in either case the 
value chosen will not maximize his satisfaction function considered 
alone. A simple analytical expression is suggested as a measure of 
freedom in this case, and some problems which arise from this sugges- 
tion are mentioned. 


The concept of freedom is one of the most important in sociology. 
If a systematically developed mathematical sociology is to be 
achieved, a precise definition and a mathematical treatment of the 
concept of freedom must be made. Inasmuch as subjective state- 
ments about more or less freedom are constantly made, a psycho- 
physical scale can certainly be established for its measurement 
(Rashevsky, 1951, hereinafter referred to as S, B.; Chapter XXIII). 
For theoretical purposes this is, however, not enough. We need a 
mathematical formulation of freedom and of its relation to other 
biosociological quantities. 
In previous publications (Rashevsky, 1940, 1947, also S. B.) we 
have outlined a possible approach to the problem, based on the 
consideration that an individual’s freedom is closely related to his 
ability to choose any activities which he likes. On this basis a 
simple relation between freedom and population density has been 
‘derived. 
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But the choice of an activity or of a set of activities by an in- 
dividual is largely determined by the amount of satisfaction which 
the individual derives from this choice. It is therefore natural to 
look for a connection between the mathematical theory of freedom 
and that of hedonistic behavior, which is determined by a satisfac- 
tion function. Various mathematical aspects of hedonistic behavior 
have been studied by N. Rashevsky (1947, 5. B.) and by Anatol 
Rapoport (1947a,b,c). This satisfaction determines not only the 
choice of an activity, but also the choice of the quantity, or in- 
tensity, of an activity. It is not infrequent that this quantity 
has an optimal value, for which the satisfaction has a maximum. 
Sometimes the existence of a maximum of the satisfaction func- 
tion is due to the circumstance that the satisfaction function con- 
sists of two components. One is positive and increases monotoni- 
cally with the object of satisfaction; the other is negative and de- 
creases with the effort which is necessary to obtain a certain 
amount of the object of satisfaction. Since this effort is usually a 
monotonically increasing function of the amount z of the object of 
satisfaction, the satisfaction function § may have a maximum with 
respect to z. 

A maximum of § with respect to may, however, be due to purely 
physiological reasons. Certain pleasant stimuli may have an opti- 
mum value. An individual may like to smoke four cigarettes a day, 
but would not like to smoke more. A sound may be increasing in 
its pleasantness up to a certain intensity and decreasing in pleas- 
antness beyond that intensity regardless of the effort involved in 
its production. We shall consider here the second case only, where 
the maximum of the satisfaction function is not determined by the 
efforts involved. 

Consider an individual who indulges in an activity which results 
in an object of satisfaction X, in the amount z,, and who is char- 
acterized with respect to that object of satisfaction by a satisfac- 
tion function $,(%). The individual will choose such an amount 


@, = xt which makes S,(z) a maximum. In other words, 2* is the 
root of 


dS ,(#,) 
a 20 
ra (1) 
We may say that the individual is free to choose any amount of 
X which he likes if he always can choose 2, = x*, that is, if he 
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always can make equation (1) satisfied. Limitations of this pos- 
sibility may be twofold. First they may be of a purely physical 
nature. For example the root x* of (1) may lie outside of the 
physically possible range of variation of z,. In that case either 
the lower bound 2‘ or the upper bound 2%’ of x will be chosen, de- 
pending on whether oS ,(z,)/dz, is larger for #, = 7‘, or for 2, = 27. 
The impossibility of maximizing S,(a,) will be subjectively felt by 
the individual as a limitation of his freedom of choice. With Rapo- 
port (1947a,b,c; also Rashevsky, S. B.) we shall assume that the 
strength of desire for an optimal value of «, is proportional to the 
absolute value | dS ,(«,)/da, |. That desire is then zero. In other 
words, the desire is satisfied, when (1) holds. If the smallest 
possible desire is not zero, then we may consider the magnitude of 
this smallest possible desire as a measure of the lack of freedom. 
The greater the smallest possible desire, the less the freedom. 


Put 
: dS .(z,) 
si(a) = ( s ‘| ; (2) 
vy x ,=a 
and 
Pee ae for Si (#4) > S4(27); a 
t., af for S4(24) < 84 (27). 


Then | S4(z,,,)| is the smallest possible desire of the individual 
for the optimal value of z,. We may conveniently choose a scale 
for the freedom F, which assigns to it the value 1 (complete free- 
dom) when (1) can be satisfied and which tends to zero as the 
smallest possible desire | S4(z,,,) | increases. One of the simplest 
possibilities is to put 
1 
Sa a fe 
Sociologically, however, a more interesting limitation of the 
freedom of choice may be introduced by the presence of other in- 
dividuals. Consider the case of two individuals, both producing 
the same objects of satisfaction X. Let the amount of X produced 
by the second individual be «,, and his satisfaction function be 
S,(#,). Let there be no physical limitations imposed now on the 
ranges of z, and z,. In general, if the two individuals interact 
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socially and are not completely isolated, they cannot avoid sharing 
the amounts of X produced. Thus if two individuals are in close 
proximity and one produces a sound of intensity z,, the other—the 
same sound of intensity #,, then both are affected by a sound of 
intensity «, + @,. It is in general impossible to divide the total 
intensity between the two individuals. On the other hand if the two 
individuals live in rooms separated by a wall, which transmits 
only a fraction p of the sound intensity, then the first individual 
will be affected by an intensity x, + pw,, while the second will be 
affected by an intensity px, +a,. Quite generally, denoting by 
f,(@,,%,) and f,(#,,2,) two functions of z, and 2, we may con- 
sider the case in which one individual is affected by an amount 


2, = f,(@1. 2g) (5) 


of the object of satisfaction, while the other is affected by the 
amount 


8, = fal@,s@,)- (6) 
While 
f1(@,,0)= 2,3  — f,(0, 2) = ay; (7) 
we have in general 
f,(0,@,) # a, 5 fa(w1,90) # ay. (8) 


For any given values of w, and @,, the satisfaction of the first 
individual is now equal to $,(2,), while the satisfaction of the 
second is equal to S,(2,). 

If both individuals act egoistically, then each of them tries to 
maximize his own satisfaction function, regardless of the other. 
In other words the first individual will choose such a value of 2, 
for which 


dS ,(2,) a dS (24) Oz, 


Oa, da, Oe mis (9) 


while the second will choose an #, such that 


OS (24) ‘3 d8,(2,) 02 


= is 
da, da, 08» Cu) 
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Let the roots of the system of equations (9) and (10) be w**and 
ws*. Let 


ai* = f,(et*, 23") | Ge 


83 = f,(ot*, 22°). 


The value S,(2**) is the maximum value of the satisfaction S, 
of the first individual when the second individual is present; 
whereas S,(23*) is the maximum value of S, when the first in- 
dividual is present. Except for a very special choice of the dif- 
ferent functions involved, we have in general 


epee; et zey, (12) 
where 2* is the root of (1), and x is the root of 
Si(e,) =0. (18) 
From (12) it follows that in general 
S,(2¥*) <8, (e*); S,(2t*) < S,(#%). (14) 


Remembering that z2¥* is the amount of object of satisfaction 
affecting the first individual in the presence of the second, while 
z** is the amount of object of satisfaction affecting the second 
individual in the presence of the first, we see that the satisfaction 
of each individual is now less than it would have been if the other 
were not present. Each individual would like to increase his sat- 
isfaction to the value S,(#*) and S,(##) correspondingly, but this 
is impossible under the conditions of egoistic behavior. The de- 
sire of the first individual for an optimal value of 2, is | S4(2¥*) |, 
while the desire of the second for an optimal value of 2, is 
| S4(2**) |. Hence their freedoms now are 


1 1 
at ‘< —— 15 
a 1+ 842(2**)’ 2 1+ 83% 23") Sa 


Thus, without any physical limitations imposed on the ranges of 
variation of #, and 2,, we find now Fi <1; Fj <1, while for iso- 


lated individuals F, = 1; F, =1. 
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If the interaction of the two individuals is characterized by the 


matrix (S. B., Chapter XIV) 
1 ky 
, (16) 
k 1 


then in the presence of both individuals, the first individual tries 
to maximize 


8, =S, + hes? (17) 
while the second tries to maximize 
8S, =k,,8, +8,- (18) 


We thus have the system of equations 


dS, dS, 2, , dS, a2, 


Ses + —— ———S- = . 19 
dz, da, da, afl EB dz, : (19) 
Se ren 
dx, a0 da, da, G4, dae 


Let z, and z, be the solutions of (19) and (20), and let 
2, =f (#1) %Q)3 2, =fo(%1> 2)- (21) 
The satisfactions of the two individuals are now 
S(3,. 34) = 8,(2,) + &,,8a(4,)3 (22) 
S, (2s. 3g) = yp 8(2,)4 S,(2,)- (23) 


The satisfaction S ,(Z,) of the first individual taken alone is now 
again in general less than S,(x*); similarly S,(Z,) <S,(#¥). The 
freedoms of the two individuals are accordingly given by 

i 1 
oe F 


1 et + 843(z,)’ ae fs s’7(z,) re 
Both F’;’ and F;,’ are functions of the parameters &,, and k,,, and 
for k,. =k, = 1 the case reduces to that of an altruistic behavior. 

We may, however, interpret equations (19) and (20) in a different 
manner: The first individual maximizes (17) and the second—(18), 
because actually their respective satisfactions are given by S, and 
S,. In other words the satisfaction of an individual in the pres- 
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ence of another is equal to the satisfaction of the first when iso- 
lated plus a constant multiplied by the satisfaction of the other 
individual when the latter is present alone. 

With this interpretation it may very well happen that Sq CAA es 
S (zt) and S,(2,, 2,) > S,(z%). In this case the freedom of each 
individual may increase when the two interact. 

Such an interpretation of § is different from the one which we 
have adopted hitherto (S. B.). 

It would be of interest to investigate possible relations between 
Fy and Fy’; F’’ or Fj and Fj’. Could it be that always Fe Ss 
F3’ > Fy? And would k,, = k,, = 1 give the highest possible value 
of F{’ and F3’? Those and similar problems, as well as a generali- 
zation to the case of several independent activities and for the 
case of N individuals, are suggested by the present approach. 

A different type of limitation of freedom is possible due to a 
restriction of the total amount of possible objects of satisfaction. 
Thus suppose the object of satisfaction is the spending of a given 
fraction of the day in driving a car, and let the optimal fraction for 
the first individual be x*, for the second x%. If only one car for 
the two is available, then it may happen that‘a simultaneous maxi- 
mizing of the satisfaction function of both individuals is impos- 
sible. The value z* is determined by (1), the value of w} by 
S3(z@,)=0. But ,+a%,=1. If in the z,,az, plane the point 
z*,2* lies above the line x, + 2, =1, then obviously both in- 
dividuals cannot maximize their satisfaction. What will happen? 
They will agree on some point which lies on the line z, + 2, =1, 
because for any point below both could increase their satisfactions 
by moving closer to the line. But on what point of the line will 
they agree? The simplest assumption to make is that they agree 
on such a point x4,a4 at which | S4(2’) | = | $4(@%) , that is, when 
their desires for a maximum satisfaction, and therefore their free- 
doms, are equal. More generally we may assign to the two in- 
dividuals ‘‘personality factors’ p, and p,, and assume the equal- 
ity of the products 


p, | 84(24)| = p_185(#5)1- (25) 


The factors p, and p, may represent pure physical strength. 
The actual willingness of an individual to obtain something by 
force will depend upon his strength as well as upon his desire for 
obtaining his goal. An equilibrium may be assumed when the prod- 
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ucts of the two are equal, and this leads to (25). But p, and p, 
may also represent the ability for persuasion without the use of 
physical force. 

The freedoms of the two individuals are now 


1 1 
= ——_______} = —_,—.. (26) 
TUNES 8?(«4) 2 1 + 837(#3) 
Because of (25), we may write 
Z 
F, = ET Se . (27) 
1+ 1 847(24) 


2 


If p, >p,, then F, >F,; if p, <p,, then F,<F,. Thus the 
physically stronger or the more persuasive individual has in this 
case a greater freedom. 

The value of 2 is determined by the equation 


84(0;) = S4(1-2%), (28) 
and wv, follows from #4 + 2) = 1. 


This work was aided by a grant from the Dr. Wallace C. and 
Clara A. Abbott Memorial Fund of The University of Chicago. 
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